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ABSTRACT 



A comprehensive study of multiresolution signal processing techniques is conducted. 
Background material in functional analysis and Quadrature Mirror Filter (QMF) banks is 
presented. The development of Mallat’s algorithm for multiresolution decomposition and 
recon structi on is outlined and demonstrated to be equivalent to QMF banks. The Laplacian 
pyramid and the a trous algorithm are described and demonstrated. General multiresolution 
structures are constructed from cascades of QMF and pseudo-QMF banks and are demonstrated 
for applications in signal decomposition and reconstruction and signal detection and 
identification. 
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I. INTRODUCTION 



A. MOTIVATION FOR STUDY 

In recent years a novel approach to the field of digital signal processing has received 
significant attention. This new approach, the wavelet transform, has displayed the potential for 
applications in the fields of speech and image processing, in particular. Techniques based on the 
wavelet transform have been praised for, among other things, their efficiency of computation. 

Unfortunately for the electrical engineer specializing in the field of digital signal 
processing, the concepts on which the wavelet transform have been based have been somewhat 
elusive. Pioneers of the field have been specialized in such diverse fields as geophysics, 
mathematical physics, and a branch of applied mathematics known as functional analysis. Not 
completely unrelated to the familiar Fourier transform methods, wavelet transforms extend the 
concept of signal decomposition through basis function expansion to a more general and abstract 
realm. Such abstract mathematical disciplines as set theory are invoked within the context of 
embedded vector spaces and their relationships to each other. 

Included among the purposes of this present study is an attempt to bridge the gap 
between the mathematician and engineer. It is intended to provide a rudimentary enough 
understanding of some of the concepts of functional analysis to facilitate a more in depth 
understanding. Additionally, a strong tie is demonstrated between wavelet-based methods of 
decomposition and areas within the field of electrical engineering which may be more familiar. 



Finally, an additional purpose of the present study is to evaluate wavelet transform 
methods for detection and identification applications. Wavelet processing methods will be 
applied to familiar signal processing problems such as the detection of signals in noise and the 
resolving of disparate frequency components of which a signal is composed. 

B. OUTLINE OF STUDY 

Chapter II will introduce concepts from functional analysis which commonly appear in 
the literature on wavelet transforms. Not intended to provide an in depth introduction into the 
field of functional analysis, the chapter provides a few basic tools to support further study. In 
Chapter III, the basics of multirate system theory will be explored. In particular, the Quadrature 
Mirror Filter (QMF) bank for two and more channels will be demonstrated. Chapter IV will 
introduce the theory of multiresolution analysis. From the definitions introduced, Mallat's 
algorithm for the discrete wavelet transform will be introduced and shown to be equivalent to the 
QMF bank. Two earlier multiresolution decompositions— the Laplacian pyramid and the a trous 
algorithm for the wavelet transform— will then be reviewed prior to extending the concepts of 
multiresolution analysis to general structures constructed from cascaded QMF banks. In the 
Chapter V, the basis functions of the wavelet transform and filters for QMF banks will be 
considered in greater depth. The chapter will include a discussion of two-scale difference 
equations and some basic filter design methods for QMF banks. Chapter VI will demonstrate for 
detection and identification applications the application of some multiresolution structures 
introduced in Chapter VI. Multiresolution structures will be compared with the periodogram 






decomposition from the perspective of computational efficiency, robustness with respect to 
noise, and the ability to resolve proximate spectral components. 
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II. RESULTS FROM FUNCTIONAL ANALYSIS 



A. INTRODUCTION 

In essence, wavelet transformation represents a recent development in the field of 
functional analysis. Consequently, a brief study of this branch of mathematics provides an 
appropriate starting point for a study of Wavelet Transform (WT) analysis. 

In mathematical terminology, a functional is defined as [1] "a function that has a domain 
whose elements are functions, sets, or the like and that assumes [scalar] numerical values." 
Fourier transformation--an operation familiar to electrical engineers working in the field of 
signal processing— is an example of a functional analysis technique. The Fourier transform 
integral 

f(co)=s^’{f(x)}(co) =\Z^ f(x)e" J “ x dx (2.1) 

is a functional which, through projection on a set of basis functions "e' Jt0 \" maps a continuous 

function to a specific value for a given value of "co." Similarly, the Discrete Fourier Transform 
(DFT) is an example of a discrete functional. For a discrete sequence x T =[x[0] x[l] ... x[M- 1 ] ] , 
the DFT 

X=\V M x, (2.2) 

where [W M ] k m = W M k m = e' j2ltk iaM maps the vector x to a discrete set of points X[k], 

B. NORMS AND NORMED SPACES |2| 

Linear transformations generally involve the mapping of a function (or vector) from one 

vector space to another. The study of these transformations requires the capacity to quantify 
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distances within these spaces. Measures for functional spaces possess analogies in 
geometry— such as length of a vector, distance between the points defined by two vectors, and 
the scalar product of two vectors. These measures, or norms, provide the ability to transform 
each element u in a vector space V into a real number. For a linear vector space V over the real 
number field R, the norm ||u|| for every element u € V satisfies the following conditions: 

||u|| > 0, and ||u|| = 0 if and only if u = 0. (positivity) 

||ot*u|| = |a|-||u|| for a s R. (homogeneity) (2.3) 

||u + v|| < |u| + || v|| (triangular inequality). 

Examples of norms which sometimes arise within the context of functional analysis 
follow. The most common norm is the L 2 -norm. Pythagorean theorem represents a special case 
of the L 2 -norm. For a function "x(t) M defined in the closed interval C[0, T], the 
L 2 -norm— denoted by ||x|| 2 — is defined as 

IMI 2 = [Jo ! x (t)l 2 tit] ' 2 . (2.4) 

Another norm is the sup-norm or supremum norm. The supremum of a function or 
functional is defined as the least upper bound of a function. The concept of supremum is similar 
to the concept of maximum, however a function does not necessarily ever assume the value of 
the supremum. For example, consider a function f(x)=x 2 . If the domain of f(x) includes the 
closed interval [0, 1], then the extremum (and, incidentally, the supremum) of fix) occurs for 
f(x)=l. If however, the domain of fix) is restricted to the open interval (0, 1), then the function 
fix) has no extremma. At its limit, fix) approaches unity, however it can never assume that 
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value. In this case, however, f(x) does possess a supremum, since f(x) is bounded above by 
unity. 

To applying the concept of supremum to define a sup-norm, consider the vector space 
V=C[0, t 0 ], the set of all real-valued, once-differentiable, continuous functions of t in the closed 
interval [0, t 0 ]. The sup-norm of x, HxH^ is defined as the supremum of the function "x(t) M on its 
domain, or 

||x|L = sup{ |x(t)|: 0<t<t o }. (2.5) 

The Lebesgue norm represents another norm which is related to the L 2 -norm. For a real 
number p € [0, «>) n R, the Lebesgue norm for the function u(t) defined on the interval [0, T] is 
defined as 

l|u|l p = [lo lu(t)|p dt] l p < oo. (2.6) 

If a norm can be defined in a given space, then that space can be characterized by that 
norm, for instance, if the Lebesgue norm is defined for a space of interest, than that space can 
be classified as a Lebesgue space. Furthermore, the space by which a norm is characterized is 
indicated in the subscript of the norm operator symbol. In some cases, such as for Lebesgue 
spaces, the subscript indicates the metric of the space. In others, such as sup-norm spaces, 
subscripts less indicative of the norm operation appear. In general, for some arbitrary, 
unspecified space U, the norm operator for that space is denoted by 

• 

lu 

This notation will be used throughout the remainder of this section to indicate general normed 
spaces. 
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Inner product spaces are of significant interest to signal processing. An inner product 
space is defined to be a linear vector space on which an inner product can be defined. The 
concept of inner product is related to the definition of the norm of a space. For instance, the 

L 2 -inner product (u,v) 2 is evaluated as 

(u, v) 2 =J n u(x)- v*(x) dx (2.7) 

for two vectors u, v e L 2 (f2). Unless otherwise indicated, throughout the remainder of this 
paper, inner products will be assumed to be L 2 inner products and the corresponding 
distinguishing operator subscripts will be suppressed. Additionally, inner product operations 
possess the following properties: 



1. (u, v) = (v, u) 

2. (au, v)=a(u, v) 

3. (u.+i^, v)=(u,, v) + (U 2 , v) 

4. (u, u) > 0, and (u, u) = 0 if and only if u = 0. 

5 . | (u, v) | < V(u,u)(v, v) = ||u|| ||v|| 



Symmetry 

Homogeneity 

Additivity (2.8) 

Positive Definite. 

Cauchy -Schwartz inequality. 



Given an inner product space V, the concept of orthogonality is also important for 
characterizing the space. Two vectors u, v e V are said to be orthogonal if (u, v) = 0. 
Furthermore, it is possible to partition an inner product space into orthogonally complementary 
subspaces. For instance, consider the inner product space V k . If V k _, is a subset of V k (denoted 
in mathematical symbology by V k , c V k ), then the orthogonal complement W k _, of V k ., is 
defined as 
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w k-i = v k-. 1= { u 6 v k : ( u - v ) = 0 v v€ v k-il- (2.9) 

Additionally, the union of two orthgonally complementary subspaces V k _, and W k _, to obtain a 
third is denoted V k , © W k _, = V k . The concept of an orthogonally complementary subspaces 
proves critical when defining a multiresolution analysis. 

The Hilbert Space represents perhaps the most important inner product space for signal 
processing applications. An abstract Hilbert space is an inner product space which possesses the 

following characteristics [3]: 

1 . Linearity— The operations of addition and of multiplication by real or complex numbers 
are defined for its elements: 

2. The metric of an Hilbert space is deri ved from its inner product. Consequently, for any 
two elements u and v, there is an associated real or complex number. 

3. Completenes fj -If a sequence of elements {u,,} satisfies the condition ||u„ - uj 1 — >0 V 
m, n —> <*>, then there exists an element, u, such that | |u„ - u|| — » 0 V n — » «>. 

C. LINEAR OPERATORS AND TRANSFORMATIONS 121 

Many signal processing applications of interest involve linear transformations from one 

linear vector space U to another linear vector space V. If for some vector u in inner product 
space U, there is a corresponding vector v in image space V, then the operation T, by which u 
and v are related, constitutes the corresponding transformation. Mathematically, a mapping by T 
from space U to space V is denoted by T: U —> V. 
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With respect to some arbitrary linear transformation T: U V, the following properties 

apply: 

1. T{au, + (5u 2 } = aT{u,} + PT{u 2 } Linearity 

2. ||T{u,} - T{u,}|| v < M ||u, - u,|| L . for some arbitrary M>0 Continuity 

3. ||T{u}|| v < M ||u|| L , for some M > 0 Boundedness from Above (2.10) 

4. ||T{u}|| v > C Hullu for some C > 0 Boundedness from Below 

Two important consequences for linear operators result from the properties of (2.10). The first 
states that for some invertible, linear operator T: U — > V which is bounded from below, 

l|T* , {v}|| ^ l|v|| v /C. (2.11) 

Equation (2.1 1), proven in [2], expresses the Bounded Inverse Theorem. In other words, the 
inverse of a linear operation bounded from below is bounded from above. 



The second significant result stems from the linear property. Given T: U — » V, for 
finite-dimensional vector spaces U and V, constituent vectors u and v, respectively, and bases 
{<(>,, <J) 2 , ..., 0 n } for U and {y,, y : , ..., y m } for V, then u and v can be represented as 



u= I oc k <J> k 

k=l 

m 

V = XpicVk 

k=l 



( 2 . 12 ) 



Now, writing down the form of the transformation, T: U — » V, the terms in (2.12) are related by 

v = I p jVj = I a k T{(J> k }=T{u}. (2.13) 

j=l k=l 

Since € U and € V, the mapping of {(J>j — » {y^ can be expressed as 
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(2.14) 



m 

T{^k } ~ X tj k \|/j . 

J=1 

In vector-matrix notation, (aj, {fy} and (tj k } are related by 



p> 




* 1,1 * 1,2 • 


* 1 ,n 




a , 


p. 




* 2.1 t 2,2 • 


* • * 2.11 




<*2 


_ Pm 




*m.l *m,2 


*m,n 




cx n 



Equation (2.15) represents a linear transformation of the projections onto the bases of U and V. 

In signal processing, and other sciences involving representation in terms of linear 
transformations, mappings which are one-to-one are of significant interest. Given T: U V, if 
for each u e U there exists a unique v e V, then the transformation T is defined to be 
isomorphic. Furthermore, for T{u} = v, given an isomorphic operator T, there exists a unique 
inverse operator T 1 such that u=T' {v} [4], Isormorphic operators are particularly useful 
because they can be inverted. With respect to signal processing, if a signal is decomposed by an 
invertible operator, then it can be reconstructed to recover the original signal or process. 

Before the final concept can be defined, elaboration on the definition of a functional is 
necessary. Linear functionals consist of a subset of linear transformations T: U -» V in which 
the transform space V represents a scalar field [5]. In the case of functionals, the transform 
space V represents the dual or conjugate space of U. The dual space of U is assigned a special 
notation U*. The notion of a dual space differs from that of a biorthogonal, or Riesz space 
which shall be addressed in a later section. 
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Finally, the concept of adjoint operators occasionally plays a role in functional 
representation of sequences. Given a bounded, linear operator T: U — » V for normed, linear 
spaces U and V, a linear functional § can be defined such that <j>( v )=(|>(T{u})=£(u). The resultant 
functional (--for u— is linear since <\> and T are linear. It is possible, therefore, to define an 
adjoint operator T\ such that ((u)=((T*{v })=<}>( v). The functional (, by definition, lies in U\ 
Additionally, omitting the argument vectors, the transform T* can be expressed as (=T‘<1>. In the 
case of Hilbert spaces, an operator and its adjoint are related by inner products. Specifically, for 
u, v e X a Hilbert space, 

(T{u(,v) = (u,r{v)). (2.16) 

D. REPRESENTATION JN INNER PRODUCT SPACES 

Representation theory consists of the theory of representing sequences or sets in terms of 

projections upon sets of vectors. The Riesz Representation Theorem provides the foundation for 
representation theory for functionals. If i represents a bounded linear operator in an Hilbert 

space X the Riesz representation theorem [2] states that there exists a unique vector v 0 e K such 
that 

(( w) = (v 0 , vv) V W € X 
This vector v 0 is called the representation for operator (. 

The Fourier Series Theorem [2] provides concepts fundamental to representation of a function in 
a Hilbert space. For a countably infinite, orthonormal vector set {un}”^, e X then a series of the 



form £ ^ D u n converges if and only if £ I^J 2 < 00 . In this case, the series £ snU n converges 

n=0 n=0 n=0 

to the same element x irregardless of the ordering of the terms. Furthermore, the element x, the 

orthogonal set {u n } and the weighting factors {£„} are related by 

^„ = (x, u n ) (2.17a) 

and 

x = £ £ n u„. (2.17b) 

n=0 

By linearity, Parsevals Equality follows directly from (2.17): 

(y, xj = £ (x, u n )^y, u n j. (2.18) 

Finally, if x = y, another relationship is obvious: 

I x ||‘ =(x, x)= £ |(x, u„)| 2 . (2.19) 

At this point, consideration of additional concepts is justified. The concept of closure 

arises in representation theory [2], A subset S of a normed space is said to be closed if it 

contains all its limit points. If S is a closed set, and if a sequence converges, then the limit to 

which the sequences converges is contained in S. The concept of denseness is related to closure. 

If S is an open set (not closed), the set of all additional points necessary to obtain a closed set 
including S is denoted by S. Consequently, the union of the sets SuS results in a closed set. A 

space S is said to be dense if for any vector v in S, there is an element in S which is arbitrarily 

close to v. Practically speaking, within a dense space S, it is, therefore, possible to define a 

representation with arbitrary precision if enough basis vectors are employed. 

Finally, in representation theory, the concept of completeness of an orthonormal set 
represents an important concept. Orthonormality, in the usual sense, refers to a set of vectors 
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{u n }n=i such that (u n , uj = 5 mjl where 6 mjl represents Kronecker's delta function. If an 
orthonormal basis {u n }n=i € % is complete, then the only vector v € K such that (v, u n ) = 0 V 

integer n € [1, N] is a vector such that ||v|| = 0. 

In representation by summing projections on a set of orthonormal basis functions U, one 
of three cases can exist [6]: 

1 . The set of basis functions is incomplete. 

2. The set of basis functions is complete. 

3. Some subset U s of the basis functions constitutes a complete set where the 
complementary subset U s * 0. 

In case 1, it is not possible to obtain a complete representation from U. In other words, 
(2.17) does not hold. Since U is incomplete; for the set of basis functions {u n }„ =1 , there is 

some non-zero vector v such that (v, u n ) = 0 for any integer n on [1, N], The set of Rademacher 

functions [7] represents one such example. 

In case 2, a complete representation is obtainable and (2. 17) does hold. Classical Fourier 
expansion represents an example for case 2. 

The third case is less common but does arise occasionally. In the third case, some 

arbitrary basis vector u n may lie in the closed linear span of all others in the set: 

u n = I a k u k . (2.20) 

k*n 

Such an occurrence is referred to as a frame. In the case of a frame, (2.17) does not hold. Due 
to (2.20), evaluation of (2.17b) would result in a representation of vector x containing redundant 
information. Representation by Fourier expansion with overlapping subdomain basis functions 



13 



exemplifies the use of a frame. Specifically, representation by expansion on integer translates of 
polynomial splines is employed for some signal processing applications [9]. 

The concept of a frame originates from a theorem outlined by Riesz and Sz.-Nagy 
addressing biorthogonal systems [3]. The present term, however, was introduced by Duffin and 
Schaeffer in work related to nonharmonic Fourier analysis [8], In contrast to expansion by sets 
of complete orthonormal bases, expansions by frames do not converge to a specific vector. 
Instead, such expansions converge to a specified range: 

A ||x|| 2 < Z I (4>k, X )| 2 <B||x|| 2 V k € Z (2.21) 

k 

where Z represents the set of all integers. The constant factors A, B € R are called the frame 
bounds. In order the operation to be invertible (or even, possibly, nontrivial), the requirement 
exists that A > 0 and B < °°. Daubechies [6] presents detailed descriptions of the procedures for 
calculating frame bounds for various situations. 

Employing concepts from linear operator theory, the projection operation, T, of some 
vector, x onto a set of basis functions V m e Z is characterized as a mapping from a 
Hilbert space % to the sequence of all square summable sequences l 2 ( Z): 

T: X-> ( : (Z). 

Specifically, the projection operation T, defined to be the frame operator, is defined as 

(T{x}) k = (<{> k , x). (2.22a) 

The adjoint frame operator, T\ results in the expansion on the basis set {<J) k } 

T*{c} = Xck <{>k • (2.22b) 

k 
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Combining (2.22a) and (2.22b), (2.21) can be expressed in its most abstract form as 

Al < T’T < Bl 

where 1 denotes the identity operator. 



A decomposition by frames can be reconstructed exactly if a biorthogonal system , also 
occasionally referred to as a Riesz basis [9], is constructed. The form of such a reconstruction 
appears as 

x = I(<j> k) x)<j> k =X(4> k , x)(}) k . (2.23) 

k k 

In (2.23), the vector set {<j> k } is referred to the biorthogonal basis of {<J> k }. The basis function for 

the biorthogonal set is related to the fundamental set by 

<t>k = (T* T) -1 <j) k (2.24) 

Occasionally, in literature, the Riesz basis will be referred to as the dual frame. This term 

should not be confused with the concept of dual or conjugate space previously discussed. 



Daubechies [6] describes and proves three generalities regarding biorthogonal systems. These 
include: 

1. The family (<j) k }v k with (j> k = (T*T) -1 <{> k constitutes a frame with bounds B' 1 and A' 1 . 



2 . 



The frame operator T associated with (J> k is given by T = T(T* T) 1 and satisfies 

T*T = (T*T)-‘ 

and (2.25) 

f * T = T* T=l, 



or, equivalently 



(T*)' 1 = T. 



3. And finally, 



T (T* jy l * jj* 



(2.26) 
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III. A SURVEY OF MULTIRATE SYSTEM THEORY 



A. INTRODUCTION 

Multiresolution analysis implementations frequently reduce to structures composed of 
multirate system building blocks. Consequently, prior to addressing the construction of 
multiresolution structures, multirate system theory will be briefly reviewed. Section B will 
address basic multirate system operations and representations. The simplest multirate 
svstem--the two-channel Quadrature Mirror Filter (QMF) bank— will be discussed in Section C. 
In Section D, the results from the preceding sections will be extended to an arbitrary number of 
channels. 

B. BASIC MULTIRATE SYSTEM OPERATIONS 

Multirate systems are comprised of three fundamental operators [10]— decimators, 

expanders, and linear (usually Finite Impulse Response (FIR)), digital filters. In this section, 
time- and Fourier-domain consequences of decimation and expansion will be demonstrated. 

Decimation consists of subsampling a discrete sequence, retaining only samples at 
integer intervals. Figure 3.1 presents the block diagram symbol for a decimator. In 
mathematical 




► 



Figure .?. /--Block diagram for M-fold decimator. 
notation, decimation of a discrete sequence by a factor of "M" is denoted by 



x iM (n)=x(M-n). 
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To illustrate the consequences of decimation, consider the bandpass, transient sequence 
plotted below in Figure 3.2. Figure 3.2 displays a lineplot of a 128-point sequence constructed 
from the real part of an exponentially-modulated, Kaiser windowed, sine (sine-over-argument) 
pulse: 



Re{s(n)} = 



sin ( rcn/8 j 






Re{e jnn/3 }. 






(3.2) 




Figure 3. 2-Time plot of discrete sequence described by equation (3,2). 




Figure i.i-Plot of sequence obtained from decimation by a factor of two of sequence (3.2). 
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Decimation by a factor of two, s i2 (n), of the sequence in Figure 3.2, produces sequence 
plotted in Figure 3.3. The length of sequence s i2 (n) is one half the length of sequence s(n). 
Furthermore, the sequences are related by 



s i2 (n)=s(2-n). 



(3.3) 



To consider the Fourier-domain consequences of the decimation operation, the Z-transform of 
the decimated sequence must first be evaluated. For the decimated sequence, 



Sx 2 (z) = £s(2n)z " 

n 

= I s(n) z -0 ' 2 



(3.4) 



Evaluation of the second summation in (3.4a) requires the definition of a sequence Ti : (n) whose 
values are unity for even elements and zero for odd elements [11]: 

ri 2 (n) = [l+(-l)"]/2 (3.5) 

Inserting (3.5) into (3.4) produces 

Si-dz) = I s(n)z -n/2 

even n 

= yZs(n) (l +(-l) n lz -n/2 . (3.6) 

= |[S(z'' J ) + S<-z'' ! )] 

The effect of decimation on the spectral content of the sequence s(n) is obtained from evaluation 
of S i2 (z) on the unit circle: 

S± 2 (e j “) = y [ S(ej w/2 ) + S(-e j <D/2 )] 

= 4-[S(e JW/2 ) + S(ej (0W ' 2 rt) '' 2 )] 

Decimation, therefore, causes the frequency spectrum of a frequency to become spread become 
by a factor of two. Additionally, the location of a spectral peak or any other distinguishable 
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feature will be translated by a factor of two with respect to its original spectral location. The 
Fourier transform indicated by (3.7) is a 4'Tt-periodic member of L 2 ([0, 4-^]). 

1 — Original bandpass sequence 

— Decimated bandpass sequence 

\ l \ 

1 l i \ 

1 l l \ 

' '' 

' ' ' 

A l* i * L_j ii 1 , 

0.5 1 1.5 2 2.5 3 3.5 4 

Oigital frequency (multiple of tt) 

Figure 5.4-Plots of magnitudes of S(e j0> ) and 512 ( 6 *“). 

Figure 3.4 compares the Fourier transforms of a sequence S(e*“) and its decimated 
version S i2 (e i “). The plots in Figure 3.4 illustrate the test sequence power distribution relative 
to its sampling frequency. The changes to sequence power distribution which occur from 
decimation are a consequence of a change in the equivalent sampling frequency. In discarding 
each separate sample, as is done in the decimation operation, the resultant sequence is equivalent 

to sampling the original sequence at a sample frequency of f si , = 0.5-f $ . 

As predicted by (3.7), decimation causes a spreading of a bandpass about its center 
frequency. Furthermore, distance between the center frequency of the decimated sequence 
and co=0 has been doubled with respect to that of the original sequence. In the example 
illustrated, the original sequence s(n) was constructed to be an analytic, bandpass process whose 
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center frequency f c was related to its sampling frequency f s by f s = f/6 and whose bandwidth 
B=0.125-f s . The decimation operation shifts the center frequency = f/3 and enlarges the 
bandwidth to B i2 = .25 f s or by a factor of two. The spreading and shifting effects occur because 
of the subsampling nature of the decimation operation: Decimation results in a change in the 
sampling frequency. The sequence s i2 (n) is equivalent to s(n) sampled at one half of its original 
sampling frequency. The peak magnitude of S i2 (e’“) is one half the magnitude of S(e i “) because 

of the factor of one half introduced by the sifting sequence (3.5). 



*• 




*• 



Figure 5. J— Block diagram of an expansion operator. 




Figure 3.6— Discrete plot of expanded version of s^,(n). 

The operation of expansion (or upsampling) represents, in a sense, the inversion of the 

decimation operation. Figure 3.5 presents the signal processing block diagram symbol for an 
expander. The mathematical notation for a two-fold expansion operation is given by: 
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st 2 (n) = 



(3.8) 



| s(n/2) V even ’n’ 

{ 0 V odd V 

= +(-l) n j s(n/2) 

Consequently, expansion entails the insertion of a zero between each sample in a sequence. 
Figure 3.6 presents a discrete line plot of the expansion of the sequence s i2 (n). Expansion 
restores the length of the decimated sequence to that of the original sequence. However, all 
odd-number elements of the sequence s t2 (n) are zero. 

Considering the Z-transform and Fourier transform of the expanded sequence provides 
further insight into the relationship between the expansion and decimation operations. The 
Z-transform of s^n) is evaluated as: 

St2(z) = 



Consequently, in the Z-domain, expansion of the decimated sequence s i2 (n) produces 

[S i2 ] T2 (z) = j[S(z) + S(-z)]. (3.10a) 

In other words, expansion of a decimated sequence reproduces the original sequence plus an 
additional term. This term is referred to as the aliasing term [12], Evaluating (3.10a) on the 
unit circle yields: 



ji(i + (-ir)s(i)z-° 

^l(l+(-l) 2n js(n)z- 2 " 



Is(n)z 



S z 



-2 n 



(3.9) 
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(3.10b) 



[Si 2 ]T 2 (e J “) = y[S(e J “) + S(-e J “)] 

= ^S[(e»®) + S(ej ( “- R >)] 




Figure 3. 7-Magnitude plots of S| 2 (e* °) and [S| 2 ]f 2 (e J w ). 

The aliasing term, therefore, consists of the original term shifted by n radians. The magnitude 

plots of S i2 (e i “) and [S i2 ] T2 (e i “) presented in Figure 3.7 illustrate this occurrence. The left side 

of the plot illustrates the restoration of the decimated sequence to its original center frequency 
and bandwidth while the right side of the plot demonstrates the generation of an aliasing term. 
The aliasing term illustrated in Figure 3.7 can be reduced through the use of a linear, lowpass 
filter. Such a filter is often referred to as an intei-polation filter. 

Figure 3.8 illustrates the time-domain results of applying an interpolation filter to the 
sequence plotted in Figure 3.6. The sequence in Figure 3.8 differs from the original sequence by 
a normalized mean-squared error of -1 1.52 dB where the normalized mean-square error is 
defined as 
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(3.11) 



N-J ^ \ 2 

E l«»l ! I.^n, 

n=0 

The interpolation filter whose response is shown in Figure 3.9 was constructed using a procedure 
outlined by Vaidyanathan [12]. A zero-phase, half-band, FIR filter was constructed using the 
McClellan-Parks algorithm and its characteristic polynomial factored. The minimum-phase 
zeros of the original polynomial were expanded to form an new characteristic polynomial for the 
filter whose frequency response appears in Figure 3.9. 




Figure 5.5— Superimposed plots of original sequence Re { s(n) } and of real pan of interpolated, expanded, 

decimated sequence, Relho’tSjJtiO 1 )} where k^ii) represents the impulse response of the interpolation 

filter. 



The source of the error between the original sequence and its interpolated versions is 
readily apparent in Figure 3.9. The spectral content of the interpolated sequence almost exactly 
coincides with that of the original sequence in the region below the Nyquist frequency. 
However, since the interpolation filter consisted strictly of real coefficients, its frequency 
response included an image in the half of the spectrum above the Nyquist frequency. The 
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location of the image coincided with the location of the aliasing term generated by expansion. 
As a result of partial transmission of the aliasing component, the interpolated sequence whose 
spectrum is plotted in Figure 3.9 is no longer an analytical sequence. This accounts for the loss 
of symmetry evident in the plot of the interpolated sequence presented in Figure 3.8. 




Figure 3. 9— Superimposed plots of spectral content of original sequence and interpolated sequence and of 

interpolator frequency response. 

C. TWO-CHANNEL QUADRATURE MIRROR FILTER BANKS 

The primary objective of the present section is to introduce the Quadrature Mirror Filter 

(QMF) bank. QMF bank theory is a well-developed area within the field of electrical 
engineering and can be understood with a grasp of basic linear systems theory. The purpose for 
the study of QMF banks in this chapter is to prepare for the construction of multiresolution 
structures in the next chapter. In Chapter IV, Mallat's algorithm for the discrete wavelet 
transform will be shown to be exactly equivalent to the QMF bank. Thereafter, multiresolution 
analysis structures will be constructed from QMF banks of various numbers of channels. 



x ( n ) 

► 




x ( (n) Vj (n) u , ( n ) 

Figure 3.10 - Block diagram of two-channel, Quadrature Mirror Filter (QMF) bank. (After [12]) 

Having considered the basic building blocks of a linear, multireate system, attention in 

this section will be directed towards assembling those components into a basic system. 

Two-channel, quadrature mirror filter (QMF) banks represent the most basic structures for 

transmultiplexers, sub-band coders, and discrete wavelet transforms. Figure 3.10 presents a 

block diagram of a basic, two-channel QMF bank. The vertical, dotted line in Figure 3.10 

divides the system into analysis and synthesis banks. Each of the filters in the structure is a 

half-band filter. The analysis section divides the signal, by frequency, into two channels. The 

synthesis bank then recombines the channels and generates an approximation x(n) of the original 
signal. In general, x(n) will differ from the original sequence because of three sources of error 

[13]: aliasing, amplitude distortion, and phase distortion. The Z-transform analysis presented in 

the previous section permits precise characterization of these errors. 



If a filter bank is designed such that each of these errors is exactly cancelled, then the 
structure is called a perfect reconstruction filter bank. More specifically, the perfect 
reconstruction criterion is expressed mathematically as 

x(n) = c • x(n-no) (3.12) 
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for some non-zero constant c and positive integer n 0 . In other words, for a perfect reconstruction 
filter bank, the reconstructed sequence will differ from the original sequence only by a constant 
factor and a delay. 



Beginning Z-transform analysis of the structure illustrated in Figure 3.1 1 [12], the output 
of the synthesis filter for channel k (k = 0, 1) is expressed as 

X k (z)=H k (z)X(z). (3.13) 

By (3.6), the decimator output is described by 



V k (z) = j[X k (z 1/2 ) + X k (-z 1/2 )] 

= y[H k (z 1/2 )X(z 1/2 ) + H| C (-z 1/2 )X(-z _l/2 )] ' 



(3.14) 



Applying (3.9) and (3.10), 



Uk(z) = V k (z 2 ) 

= j[H k (z)X(z) + H k (-z)X(-z)] 



(3.15) 



describes the decimator output. Finally, the reconstruction component for channel k is given by 

(3.16) 



X k (z) = F k (z) U k (z) 

= y F k (z)[H k (z) X(z) + H k (-z) X(-z)] ' 



Since the reconstructed signal is simply the sum of the synthesis filter bank channels, the overall 
transfer function for the filter bank, in matrix form, becomes 

*»-« X< -) ][ h/S I p!S 1 

Now, (3.17) can be expanded and expressed as 

X(z) = \ ( F 0 (z) H 0 (z) + F , (z) H , (z ) ) X(z) + ± ( F 0 (z) H 0 (-z) + F , (z) H , (-z) ) X(-z). 
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When written in this manner, two categories of components are evident: the desired 
reconstructed components composed of X(z) and its factor, and an aliasing term consisting of 
X(-z) and its factor. Evaluated on the unit circle, the spectrum of the aliasing terms consist of a 
replica of the spectrum of the original sequence only shifted by n. 



In order to separately represent the components in the QMF bank output in terms of the 



reconstructed signal component and an undesired aliasing component, it is possible to represent 
(3.17) as a double-input, single output transfer function: 



’ T(z) ' 


_ J_ 


H 0 (z) 


H,(z) 


' Fo(z) ' 


_ J_ 


F 0 (z) H 0 (z) + F i (z) H j (z) 


. A < z ) . 


2 


_ Ho(-z) 


Hi(-z) _ 


. F '( z ) . 


” 2 


F 0 (z) H 0 (-z) + F i (z) H i (-z) 



To satisfy the perfect reconstruction criteria, it is necessary that 

A(e jw )=0 Vo (3.19a) 

and 

T(e j(0 ) = e-j (0n »| T(e jw ) |=c-e- jlDn ° (3.19b) 

where c constitutes a positive constant and n 0 is a positive integer. 

N 

Vaidyanathan [12] points out that, given a filter Ho(z) = Z ho(n)z n , the remaining 

n=0 

filters required to construct a two-channel, perfect reconstruction QMF bank are obtained from 



H,(z) = -z- N Ho(-l/z*) 

F 0 (z) = z- N H*(l/z*) 

and 

F,(z) = z- n 'h;o/z*) 



(3.20a) 



coupled with the constraint that 



H 0 (ej w ) | 2 + | Ho(ej t< ^ ) ) 




(3.20b) 
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for some positive constant c (not-necessarily the same constant as in (3.19b)) satisfy the perfect 
reconstruction criterion. A set of filters characterized by (3.20b) is said to be power 
complementary. Given (3.20), the problem of designing a two-channel, perfect reconstruction 
QMF filter bank reduces to the problem of selecting a valid analysis filter h 0 (n). Some of the 
details of selecting good filters for these applications will be addressed in Chapter V. 



At this juncture, it is appropriate to introduce some additional concepts. The matrix 



H(z) = 



H 0 (z) H,(z) 

Ho(-z) H i (— z) 



(3.21) 



constitutes the alias compensation matrix. If 

H H (ej“) H(e jco ) = dl, (3.22) 

where the superscript H denotes the transpose-conjugate of a matrix, I represents an identity 
matrix, and d is a positive constant, then the matrix H(e’“) is said to be unitaiy. By definition, 
unitary operators are operators which, when applied to their inverses, produce an identity 
operation. Within the context of linear algebra, a matrix is unitary if each of its columns is 
linearly independent from all of the other columns in the matrix. Consequently, for a QMF bank 
with an alias compensation matrix which is unitary, the aliasing component of the output is 
linearly independent of the reconstructed component and the two components are, therefore, 
separable. 

— A 

Using the tilde notation introduced by Vaidyanathan [12] where H(z)=H H (l/z‘), and 
where each element of H(z) is stable, if 

H(z) H(z) = dl (3.23) 
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then H(z) is said to satisfy the paraunitary property and the system characterized by H(z) can be 
described as a paraunitary system.. A paraunitary Z-transform matrix will be unitary when 
evaluated on the unit circle. Furthermore, if the coefficients of H(z) are all real, then H(z) is 
lossless bounded real. Additionally, Vaidyanathan lists three important properties for 
paraunitary systems: 

1 . The determinant of a square, FIR, paraunitary system produces an allpass polynomial. 



2. Paraunitary systems are power complementary. If h(e i “) = [H^e 1 ”) H^e* 01 )] 1 



3. The submatrices of II(z) are paraunitary. 

Another concept which becomes useful in multirate system theory is the polyphase 
representation. This representation will be introduced for a two-channel system here and 
extended in the next section to a structure entailing an arbitrary number of channels. It is 
possible to express the Z-transform of a sequence h k (n) as 



That is, 



det{H(z)} = az K , K > 0, a * 0.. 



(3.24) 



H 0 (e i “)p + |H 1 (e , “)| 2 = h(e j(0 )h(ej 0> ) =c V to 



(3.25) 



Hk(z) = Ih k (n)z " 



D 



(3.24) 



= Zh k (2n)z 2n + z 1 Zh k (2n+l)z 2n 



n 



n 



If subsets of the sequence h k (n) are defined as 



e k .o(n) = h k (2n) 
e k .i(n) = h k (2n+ 1) 



(3.25) 
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then it is possible to represent the two-channel, multirate system described in this channel in an 
additional manner. A z-dependent column vector h(z) can defined such that 



H 0 (z) 
H i (z 



Eo.o(z 2 ) 

Ei.o(z 2 ) 



Eo.t(z 2 ) 
E i,i (z 2 ) 



1 



= E(z 2 ) 



1 

„-i 



(3.26a) 



where 




Eo.o(z 2 ) 

E,.o(z 2 ) 



Eo.i(z 2 ) 
E (z 2 ) 



Furthermore, upon comparing (3.26a) and (3.21), the following relationship arises: 

H T (z) = [ h(z) h(-z) ]. (3.26b) 

Therefore, in terms of the polyphase matrix E(z 2 ) of (3.26a), (3.26b) becomes 

1 1 



’ h(z) h(-z) ] = 



z 1 -z 1 



E(z 2 ) = 



1 0 

0 z _1 



1 1 
1 -1 



E(z 2 ). 



Consequently, a simple relationship between H(z) and E(z 2 ) exists 

H T (z) = E(z 2 ) 



1 

o 


’ll" 


_ 0 z- 1 - 


1 _1 



(3.27) 



Obviously, the right-most matrix on the right-hand side of (3.27) is simply the transformation 



YV, matrix for a two-point DFT. The middle matrix is simply a diagonal matrix of delays and 
can be expressed as 



D : (z) = 



1 0 
0 z-‘ 



Finally, it is not difficult to demonstrate that if the alias compensation matrix is 
paraunitary, then so is the polyphase matrix: 
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(3.28a) 



H(z) H(z) = E((z- 2 )*)D 2 (1/z)W 2 W 2 h D 2 (z)E t (z 2 ) 

= E*((z“ 2 )’) E T (z 2 ) 

= dl 

or equivalently, 

E(z 2 ) E(z 2 ) = d I. (3.28b) 

Equation (3.28a) follows since W 2 W, H = I and D 2 (z) D 2 (l/z) = I. Equation (3.28b) follows 
from (3.28a) by making the substitution £ = z 1 followed by evaluating the complex conjugate of 
each side of the equation. 

Vaidyanathan delineates four important properties of paraunitary filter banks satisfying 
(3.20a) [12]: 



Where N represents the order of the filters from which the structure is constructed, filter 
banks satisfying (3.20a) result in perfect reconstruction with 



x(n) = 0.5 x(n-N); 

The analysis filters are power complementary and satisfy the condition 

I H,(ej“) 1 = 1 H 0 (ei (to- ' t) ) I; 



(3.29a) 



(3.29b) 



3. The synthesis filters are also complementary and, furthermore, they also satisfy 

| F k ( e j“) | = | H k (ej‘°) |; (3.29c) 

4. and, all filters have order N=2J+1 for integer J. This condition ensures even filter 
lengths. 



To demonstrate the operation of a two-channel, QMF bank, a simple, 256-point sequence 
was constructed. The sequence was comprised of two sinusoids windowed by a complicated, 
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Figure 3. 1 1— Time-domain plot of 256-point test sequence generated using (3.30). 




Figure 3. 72-Plot of normalized power spectral density of 256-point test sequence generated using (3.30). 
exponential window: 

s(n) = ^1 -e~ n/64 j ■ e -(n-94) '- 1 ■ e 2n/75 • sin ^-^-n j + jcos ^-^-n j . (3. 

A time plot of the 256-point test sequence (3.30) is presented in Figure 3. 1 1 and the power 
spectral density of the test sequence appears in figure 3.12. The test sequence contains harmonic 
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components at digital frequencies of 7 *tt/ 20 (= 0.175-f s ) and 1 7 -tt/ 20 (= 0.425 Q. The 
components were selected such that one appears in each half of the frequency spectum. 




Figure 3.73--bnpulse responses for 57 5 "-order FIR filters used in two-channel QMF analysis bank. 
Figures 3.13, 3.14 and 3.15 characterize the filter bank to which test sequence (3.30) was 

applied. The impulse response of the analysis filters for each channel are plotted in Figures 
3.13a and 3.13b. The low-frequency channel analysis filter, H 0 (z), a 57 <h -order FIR filter, is a 
minimum-phase spectral factor of a zero-phase, lowpass, half-band filter designed using the 
McClellan-Parks technique. From H 0 (z), the low-frequency channel synthesis filter F 0 (z) and 
both high-frequency channel filters were designed using (3.20a). 

The equivalent transfer function, t(n), whose Z-transform was defined in (3.18), is 
plotted in Figure 3.14. It has been asserted in (3.12) and (3.19b) that, for perfect reconstruction, 
t(n) will be of the form 

t(n)=5(n-N), (3.31) 
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where N is the order of the filters from which the QMF bank was constructed. Superimposed 
with a dashed-line plot, is a plot of the sequence t(n)-5(n-N). The axis for the superimposed plot 
is graduated on the right-hand side of the plot. From the superimposed plot, it is evident that the 
sequence t(n) approaches the form of (3.31) with a maximum error of less than 2 X 10' 5 . 
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Figure J./4--Superimposed plots of equivalent impulse response for QMF bank, t(n) (defined as in (3. 1 8)), and of 

t(n)-5(n-N„J. 




Figure 3. /5-Superimposed plots of frequency responses (in dB) of low- and high-frequency filters for QMF bank. 

Also superimposed is a plot of |H 0 (e ,a> )| I -t-iH l (e ,£l> )| 2 . 
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a 20 <0 63 8 a :aa 120 ua t&a 



Sample sequence number, n 

Figure 3. M— Time-domain plot of low-frequency chamiel decimator output for 256-point test sequence. 




Figure 5.77-Time-domain plot of high-frequency chamiel decimator output for 256-point test sequence generated 

by (3.30). 

The frequency responses of the filters are plotted in Figure 3.15. Only two frequency 
responses are shown because of the condition described in (3.29c). Also superimposed on 
Figure 3.15 is a plot of |H 0 (e j “)p+|H 1 (e |C) )| : \ This last plot demonstrates the degree to which the 
filter bank satisfies the power complementary property, (3.20b). 



35 



Time plots of the low- and high-frequency channel decimator outputs are presented in 
Figures 3.16 and 3.17. As a consequence of the decimation-induced translation of the spectral 
peaks for each component, the plots appear to reflect harmonics whose frequencies are quite 
close together. Additionally, due to the shape of the filter, the flat segment in the trailing edge 
of the sequence plotted in Figure 3.16 converges sharply to zero. The flat segment in Figure 
3.17 is on the leading edge of the sequence. Since the synthesis filter for channel is the time 
reversal of the analysis filter, the relative locations of the flat regions will be reversed, resulting 
in a flat segment on either side. This represents the source of the delay. The difference occurs 
since, by (3.20a), the low-frequency and high-frequency channel analysis filters are modulated 
time-reversals of each other. Also worth notice, the relative magnitudes of the low-frequency 
and high-frequency channel decimator outputs have maintained the proportions established in 
(3.30). As expected, the decimator output of the low-frequency in figure 3.16 channel maintains 
a peak amplitude greater, by a factor of four, than the peak if the envelope for the 
high-frequency channel in figure 3.17. 

Figure 3.18 presents superimposed, normalized power spectral density plots for the 
decimator outputs and the original sequence. As discussed previously, decimation changes the 

sampling frequency. The sampling frequency ^si2 for a decimated sequence is related to the 

sampling frequency f s for the original frequency by 

f si2 = f / 2 ' 
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As expected from (3.7), the peak at 7rc/20 (= 0.175 f s ) has been translated to 7 tt/ 10 (= 0.35 f si2 ). 
The peak which was originally at 17rc/20 (= 0.425 f s ) was shifted to 17 ji/ 10. Due to aliasing, the 

spectral peak located at 17 jt/ 10 appears at 20 k/ 1 0 - 17ti/10 = 3k/ 1 0 (= 0.15 f s | 2 ). Additionally, 
each of the peaks has been widened. 




Figure i. /^-Superimposed plots of normalized power spectral densities of the high- and low-frequency channel 

decimator outputs and the original sequence. 




Figure 3./9~Superimposed power spectral density plots of the original 256-point test sequence generated by 
(3.30) and the high-frequency and low-frequency channel expander outputs. 
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Figure 3.19 includes spectral density plots of the original sequence and the output of each 
of the expanders. As indicated by (3.10b), the expander reintroduces the original sampling 
frequency at the expense of creation of aliasing terms. For the low-frequency channel, spectral 

peaks occur at the original location, 7 jt/ 20 (= 0.175 f s ), and at 33rc/20 which, because of aliasing, 
appears as 13rr/20 (= 0.325 fj. Similarly, for the high-frequency channel, in addition to the 
peak at its original position of 17rc/20 (= 0.425 f s ), an aliased peak is also present at 3 tt/ 20 (= 
0.075 f s ). Although, at the expander output, the width of each peak has been restored to that of 
the input, each peak is six decibels below the original since the signal energy has been divided 
between the original term and an aliasing term. 




Figure 3.20~Superimposed plots of power spectral densities of original sequence S(e l< “) and of low-frequency 
QMF channel component and high-frequency QMF channel component . 

Superimposed plots of the power spectral densities of the original 256-point test sequence, and 
of the sequences at the outputs of the low- and high-frequency channel synthesis filters are 
exhibited in Figure 3.20. For each channel, the effects of the synthesis filters are evident. The 
aliasing terms have not been completely eliminated, but the spectral peaks of the aliasing terms 
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have been reduced in magnitude by approximately 95 dB. The restored spectral peaks, still 3 dB 

below the original, track quite closely within the passbands of their respective filters. The 
time-domain of the reconstructed sequence, s(n) = So(n) + si(n), is plotted in Figure 3.21. 




Figure 3. 27-Plot of reconstructed version of 256-point test sequence applied to QMF bank input. 

In this plot, the resemblance of the reconstructed signal to the original is evident. The 

reconstructed sequence lags the original sequence by 57, the order of the filters from which the 
QMF bank was constructed. Furthermore, if the normalized, mean-squared error is defined as 



A j s(n)-stii) | 
s : (n) 



N-i f 

X (slnHsin) 

n-" 

N— 1 

Z s-(») 

n=0 



then for the sequences plotted in Figure 3.21, e dB = -168.2 dB. 

D. M-CHANNEL FILTER BANKS 



(3.32) 



Having, in the previous section, discussed and demonstrated the implementation of a 
two-channel, quadrature minor filter bank structure, it remains to extend the results to structures 
consisting of an arbitrary number of channels. With some simple modifications to relationships 
from sections B and C, M-channel filter banks can be shown to be largely analogous to QMF 
structures. Finally, it is worth observing that in literature [13], M-channel filter banks are 
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extensively referred to as M-Channel Quadrature Mirror Filter Banks. The use of the term 
quadrature represents, in this case, a misnomer. Nevertheless, the terminology has continued to 
be applied to these more complex systems. 



The structure for an M-channel filter bank is illustrated in Figure 3.22. The structure is 
entirely analogous to the two-channel case. Each channel contains a factor-of-M decimator and 
an expander which, respectively, subsamples by a factor of M and inserts M-l zeros between 
each sample. The filters are all M^-band filters with, for the ideal case, frequency responses of 



Hk(e j “) = 



1 V £ [k-M 

0 otherwise 



(3.33) 



Furthermore, the concepts of perfect reconstruction and power symmetry also apply. 
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Figure J.22--Structure of M-channel filter bank. After [1 1]. 
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Beginning Z-transform analysis, transformation of a decimated sequence requires the 
construction of a "sifting" function equivalent to that defined in (3.5). For M channels this 



function is defined as [12] 

riM(n) = ^Z l W- nk (3.34) 

1 k=0 

where Wm = e j2n/M . Employing 3.34 produces s iM (n) = s(M-n)-T| M (M-n), thereby ensuring 
that S| M (n) is zero for non-integer values of n. For any integer n which is an integer multiple of 

M, W M "=1. Otherwise, the series term in (3.34) becomes a summation around the unit circle in 
the complex plane which is evaluated as zero. Evaluating the Z-transform of the factor-of-M 
decimated sequence s(n) produces 

Si M (z) = Zs(M-n) T| m (M- n) • z~ n 

n 

M-l f \ ~ n 

- ^? s<n) l W M z ™J ■ (335) 

= sfw m z 1/M l 

■ »ir" J 



Furthermore, the Z-transform of the decimated sequence s iM (n) after expansion by a 
factor of M is evaluated as 



[SxmW 2 ) 



Evaluating 3.35 on the unit circle produces 



Isj, M (n/M)z n 

n 




Si M (e ja ) 



JL g j e j 2 rc n\jM e j co/M 



m=0 



(3.36) 



(3.37) 



Consequently, decimation spreads, by a factor of M, the power spectral density of a sequence. 
Furthermore, the M-fold decimator produces, from a sequence with a 2 - 7 t-periodic frequency 
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response, a sequence with a 2-M-jt-periodic sequence. Similarly, on the unit circle, (3.36) 



becomes 




(3.38) 



By (3.38), it is apparent that expansion of a decimated sequence produces M-l aliased terms in 
even intervals across the spectrum. Additionally, the frequency response of the sequence has 
been restored to its 2-7t-periodocity. 

Applying (3.35) and (3.36) to the M-channel structure in Figure 3.22 produces the 



Fo(z) 

. Fi(z) 

Fm-i(z) 

(3.39) 

The matrix H(z) is (3.39) is the M-channel analogy to the alias cancellation matrix. The column 
vector containing the synthesis filter characteristic polynomials can be represented as 
f(z)=[ F 0 (z) F,(z) ••• F m _,(z)]. For alias cancellation, it is necessary that 



following matrix formulation for the system output X(z) [12]: 




H 0 (z) H,(z) 

Hq(Wmz) Hi(Wmz) 



Hm-i (z) 
Hm-i(Wmz) 



H 0 (W^-'z) z) ••• H M -,(W£J -1 z) 



42 



t(z) = 



= H(z)f(z). 



(3.40) 



t(z) 

0 

0 



Furthermore, analogous to (3.24), it is possible to define a polyphase representation of 
the characteristic polynomial H k (z): 

H k (z) 



= Ih k (n)z'“ 

= Ih k (Mn)z- Mn + z-'Ih k (Mn+ 1)z- M d + ... + z -< M+1 >Ih k (Mn + M- l) z - Mn 



M-l 

= I z~ m Ih k (Mn + m)z 

m=0 n 
M-l 

= s z- m E kjn (z M ) 

m=0 



-Mn 



(3.41) 



where Ek. m (z) ek.m(n)=hk(Mn + m). From the definition of the polyphase representation, 

the polyphase representation follows that 





H 0 (z) 




E 0 .o(z M ) 


Eo.i(z m ) • 


• E 0 .m-i(z m ) 


1 


h(z) = 


Hj(z) 


= 


Ei.o(z M ) 


E,.,(z M ) • 


Ei.m-i(z m ) 


z- 1 




Hk(z) 




E M -i.o(z M ) 


E.vi-i.i (z m ) • 


•• Em-i.m-i (z m ) 


z -<M-n 






— 


E(z M )e(z) 








where 






E 0 .o(z M ) 


E 0 .i(z M ) • 


• E 0 „m-i(z m ) 






E(z M )= 


E,.o(z M ) 


E,.i (z M ) • 


Ei,m-i(z m ) 










E M -i.o(z m ) 


Em-i.i(z m ) • 


• E M -,.m-i(z m ) 




and e(z) = [l 


[ z 1 ••• 


Z^M-D] t 









,(3.42) 



Finally, the alias compensation matrix H(z) and the polyphase matrix E(z M ) are linked by 
a relationship analogous to (3.27). First, 
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H(z) = [ h(z) h(W M z) ••• h(W^'z) ] 

E(z M )e(z) E^(W m z) m j e(W M z) e((W“-‘z) m ) e(W“-‘z) 



. (3.43) 



M 



Now, since, for any integer k, ( wj^ = ( e^ 2rtlcM = e j2rtk = 1 , (3.43) becomes: 



\ M 



H T (z) = E(z M )[ e(z) e(W M z) ••• e(W^-'z) ]. 
The block matrix on the right-hand side of (3.44) is equivalent to 

e(z) e(W M z) e(W^"'z) ] = 



( 3.44) 



1 


0 • 


0 


■ w° M 


W°M • 


w° M 


0 


z-' • 


0 


W°M 


Wm • 


•• w”- 1 


0 


0 • 


.. z -(M-l) 


. < 


3 

22 - .. 


.: w r> ! 



(3.45) 



Combining (3.44) and (3.45), therefore, produces a general relationship between the alias 
cancellation and polyphase matrices: 

H(z) = W»D m (z)E t (z m ) (3.46) 

where, 



D m (z)= 



1 0 ••• 0 

0 z- 1 ••• 0 



0 0 






Obviously, from (3.46), if H(z) is paraunitary, then E(z m ) is also paraunitary. Additionally, the 
properties expressed by (3.24), (3.25), and (3.29) are equally valid for the M-channel filter bank 
systems. 



44 




Figure 3.23 - Time plot of 256-point test sequence generated using (3.47). 




Figure 3.24 - Power spectral density of 256-pomt test sequence generated using (3.47). 

To demonstrate the operation of a three-channel filter bank system conforming to the 
structure of Figure 3.22, a 256-point test sequence similar to that of (3.30) was generated: 

s(n) = (l- e' ,v64 ) e^"- 94 ^ 1 • e 2n/75 • [cos(f£ 27 • n ) + jCos(^- 55 • n) + jcos (jfc 113 n)]. (3.47) 

Equation (3.47) employs the same envelope as (3.30) applied to the sum of three harmonic 
components. Spectral peaks in (3.47) are located at digital frequencies of 2-JT-27/256 (=0.027 f s ), 
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2-71-55/256 (=0.215-f s ) and 2 -tt 1 13/256 (=0.44-f s ). Consequently, spectral peaks occur within 
each third of the frequency spectrum below the Nyquist frequency. A time plot of the test 
sequence is displayed in Figure 3.23 while Figure 3.24 presents the power spectral density of the 
test sequence. 




Figure 3.2.5-Superiinposed plots of impulse responses of filters used to implement three-channel filter bank. 




l 



Figure 3. 26— Equivalent impulse response for three-channel filter bank structure whose filter impulse responses 

are plotted in Figure 3.25. 
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The filter bank employed for the signal decomposition and reconstruction was 
implemented using 14-point third-band filters whose coefficients are tabulated in [11]. The filter 
impulse responses are plotted in Figure 3.25. Satisfying (3.34) for the M-channel case is 
significantly more difficult than for the two-channel case. Consequently, as illustrated in Figure 
3.26, the deviation of the equivalent system impulse response deviates from the ideal case of 
(3.31) by an order of magnitude more than was observed for the two-channel structure 
demonstrated in Section C. For the two-channel structure, the root-mean-square deviation from 
(3.31) was 5.16 * 10' 6 . For the three-channel system whose equivalent response is plotted in 
Figure 3.26, the root-mean-square deviation from (3.31) is 7.42 10' 5 . Furthermore, the peak 
amplitude distortion is also greater. The peak amplitude distortion is approximately 4 x 10' 3 dB 
from the power complementary case. This represents a noteworthy increase over the 1 x 10' 3 
dB peak distortion for the two-channel structure demonstrated in Section C. 




Frequency (multiple of f ) 

Figure 3. 2 7-Superimposed plots of frequency responses for filters whose impulse responses are plotted in Figure 
3.25 (right-hand axis) and amplitude distortion from three-channel filter bank constructed from filters of Figure 

(3.25) (left-hand axis). 
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Insight is obtained from considering the decimation operation from both the time and 
frequency domains. Time-domain plots for the output of each of the decimators appear in 
Figures 3.28a, b, and c. The power spectral densities of the decimator outputs are superimposed 
in a plot presented in Figure 3.29. As discussed previously, the decimation operator, in the case 
of a factor-of-three decimation, retains only one sample out of every three. Therefore, the length 
of the decimated sequences is one third of the length of the original sequence. Furthermore, the 

effective sampling frequency ^si3 for the sequence decimated by a factor of three is related to the 
sampling frequency f s by 

fsi3 = fA 

Additionally, decimation produces aliasing terms. However, these terms are of no consequence 
since the lie outside of the region [0, 2 - k ], Figure 3.29 displays the power spectral density of the 

content of each filter bank channel with respect to its- post-decimation sampling frequency 




Figure i.2<$a--Low-fi-equency channel decimator output for 256-point test sequence generated by (3.47) applied 

three-channel filter bank of Figure 3.25. 
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Figure 3.2<$6--Medium-ffequency channel decimator output for 256-point test sequence generated by (3.47) 

applied to three-channel filter bank of Figure 3.25. 




Figure 3.28c— High-frequency channel decimator output for 256-point test sequence generated by (3.47) applied to 

three-channel filter bank of Figure 3.25. 

In the case of the low-frequency channel, the analysis filter transmits only that portion of 
the power spectral density of the test sequence which lies in the region [0, tc/3]. As discussed 
during the development of (3.36) and (3.37), because of decimation, the content of the region 
[0, 7t/3] will be linearly redistributed over the region [0, 7t]. The spectral peak passed through 
the low-frequency channel is, as a result, translated from its original location at 2 -jt27/256 
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(= 0.105-f s ) to an apparent, post-decimation location of 2-71-8 1/256 (- 0.316f si3 ). Similarly, the 
spectral peak contained in the medium-frequency channel is translated from its original location 
at 2-71-55/256 (= 0.215-Q to an apparent, post-decimation location of 2-71-165/256 (= 0.645-f s | 3 ). 
This expectation is confirmed in Figure 3.29 in which the spectral peak of the test-sequence 
component contained medium-frequency channel appears at the predicted location with an image 

appearing at 0.355- ^3- Finally, the test-sequence component passed through the high-frequency 
channel, originally appearing at a location of 2-7t-l 13/256 (- 0.414-Q, after decimation, assumes 
an apparent position of 2-71-339/256 (= 1.324-f s | 3 ). Because of aliasing, this component, in 
Figure 3.29, is indicated at 2 -tt339/256 - 2-tt256/256 = 2-tt 83/256 (= 0.324-f si3 ). 




Figure 3.29 - Superimposed plots of the power spectral densities of the original sequence and the outputs of each 
decimator for the 256-point test sequence of (3.47) applied to the three-channel filter bank of Figure 3.25. 

In Figures 3.28a, b and c, the zero-crossings of the sequences appear to occur at similar 
frequencies. This observation is confirmed by plots of the decimator output spectral densities 
superimposed in Figure 3.29. In fact, for the case under consideration, the spectral peaks of the 
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decimator outputs are separated by a maximum of 0.5-f s . This example provides insight into the 
nature of the decimation operator. In general, decimation transforms a narrowband process such 
that the result "fills" the spectrum below the Nyquist frequency of the channel. 




Figure JJO-Superimposed plots of the power spectral densities of the original sequence and the expander outputs 
for each channel of the three-channel filter bank of Figure 3.25. 



The power spectral densities for the expander outputs for each channel of the filter bank 
represented by Figure 3.25 are plotted in Figure 3.30 above. As predicted by (3.38), application 
of the expansion operator to a decimated sequence imposes two consequences. First, the 
effective sampling frequency of a decimated sequence is related to the sampling frequency of the 
original sequence by 

IWtj - 3 f si , - f, 

in the case of expansion by a factor of three. Secondly, the aliasing terms generated by 
decimation, which are originally outside of the region [0, 2-7t], are translated to within the region 



[0, 2-n], 
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In the case of the test sequence generated by (3.74), two aliasing components in addition 
to the desired component have appeared in each channel. For each channel, the aliasing 
components appear at integer translates of 2-Jt/3 with respect to the restored component of the 
original sequence. Because the spectral peaks of the original signal were separated from each 
other by roughly ji/3, each aliasing term generated for each component coincides fairly closely 
with one of the other components. This occurrence is reflected in Figure 3.30. 

In the low-frequency channel, the location of the peak of the spectral component has 
been restored to its original location of 2 *tc- 27/256=2-7C-8 1/768. However it is accompanied by 

aliasing terms at 2-tr-337/768 (=0.438-f s ) and 2-71-593/768 (=0.772-^ whose image appears at 
2 -jt - 2-TT-593/768 = 2-7t- 175/768. The restored medium-ffequency-channel component, which 
reappears at 2*tc-55/256=2-7T- 1 65/768, is accompanied by aliasing terms at 2-7T-42 1/768 
(“0.548-fj) whose image is present at 2-tt-347/768 (=0.452-f s ) and at 2-JT-677/768 (“0.882-Q 
which has an image at 2-^-91/768 (=0.1 185-f s ). Finally, to the restored component for the 
high-frequency channel which is located at 2-tt 1 13/256=2-^-339/768, are added aliasing terms at 
2-71-595/768 (=0.775-f s ) for which an image appears at 2-tt 173/768 (=0.225-f s ) and at 2-tc-83/768 
(=0.1098-f s ). 

Figures 3.3 la superimposes plots of the outputs of the synthesis filters for each channel 
of the filter bank of Figure 3.25. Each plot indicates the spectral content of the corresponding 
channel. The results of the recombination of the channels are plotted in Figure 3.31b. In Figure 
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3.31a, it is apparent that the aliasing components of the expander outputs have not been 
completely blocked by the synthesis filters. In the worst case, for the spectral region containing 
the medium-frequency channel, residual energy from an aliasing component is only 
approximately 20 dB below the desired spectral peak for that channel. However, upon 
examination of Figure 3.31b, it becomes apparent that alias cancellation does occur. The power 
spectral density of the reconstructed sequence very closely coincides with that of the original 



sequence. 




Figure 3. J/a--Superimposed plots of the power spectral densities of the original sequence and of the expander 
outputs for each channel of the three-channel filter bank of Figure 3.25. 



Figure 3.32 presents a time plot of the reconstructed sequence. Again, the reconstructed 
sequence appears to be an approximate delay of the original sequence. In fact, when the signals 
are synchronized, the normalized mean-square error (3.32) of the reconstructed signal is -66.94 
dB. The the reconstruction error exceeds that of the two-channel demonstration of Section C 
because of the poorer quality of the filters with which the three-channel filter bank of Section D 
has been implemented. As indicated by Figure 3.26, the equivalent impulse response of the 
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three-channel filter bank deviates from a pure delay by a margin three orders of magnitude 
greater than the deviation for the two-channel filter bank. 




Figure 3.3 1 ^--Superimposed plots of power spectral densities of original sequence generated by (3.47) and of 
sequence reconstructed by thee-channel filter bank of Figure 3.25. 




Figure 3. 32-Time plot of reconstructed version of 256-point test sequence applied to three-channel filter bank of 

Figure 3.25. 
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IV. THE THEORY OF MULTIRESOLUTION SIGNAL PROCESSING 

A. INTRODUCTION 

Having, in the preceding two chapters, laid the necessary groundwork, the theory of 
multiresolution signal processing will next be considered. In this chapter, Section B, presents 
the concept of multiresolution analysis. In Section C, Mallat's multiresolution algorithm will be 
developed from a projection operator perspective and the equivalence of multiresolution 
mathematical operations and two-channel QMF banks will be demonstrated. Section D will 
outline the development of the Laplacian pyramid and the A Trous algorithm, two of the earliest 
multiresolution decomposition techniques. In Section E, multiresolution structures comprised of 
cascades of filter banks will be constructed and demonstrated. 

Signal processing techniques commonly entail decomposing a signal by representing it in 
terms of its projection on a vector space. The most common method, the Fourier Transform, 
defined in Chapter II: 

f(co) = <^{f(t)}=Jl f(t)e _J0)t dt . (4.1) 

In analyzing a time-varying signal, however, (4. 1 ) presents an obvious disadvantage: Only one 
representation yectorjsjused for alljime. Consequently, time- varying aspects of f(t) are 
averaged over all time and lost. To address this shortcoming, the concept of the Short-Time 
Fourier Transform (STFT) was developed [14]: 

F((o,i)=j^ f(t) w(t-T)e" jW1 dt. (4.2) 
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Equation (4.2), through the introduction of a window function w(t), improves on (4.1). 
Typically, the window a function is either strictly time-limited or possesses a rate of decay such 
that its value outside of a limited, contiguous region is negligible. Additionally, equation (4.2) 
can be interpreted as the projection of the signal of interest on modulations and translations of a 
vector w(t). Employing concepts from Chapter II, (4.2) can be interpreted as a mapping from a 
one-dimensional space of real numbers to a two-dimensional space of real numbers, or, 

F: R R 2 . 

Additionally, the representation vector r(co, i) for the projection operator F, the operation of 
(4.2) is 

r(co, x) = w(t - x) e'-”" 1 (4.3) 

In the branch of mathematics known as group theory, the operation (4.2) belongs to a 
particular class of operators known is the Weyl-Heisenberg Group [15]. A group is a set of 
transformations satisfying the properties of closure, associativity, identity and invertibility [16, 
17], The Weyl-Heisenberg Group consists of a family of transform operators characterized by 
modulation and translation of a single representation vector [18]. 

Although it is commonly employed, two shortcomings of the STFT have been asserted. 
The representation vectors for (4.2) represent a frame in the sense described in Chapter II [6]. 
Therefore, if it is necessary invert the transform in order to reconstruct a signal from its STFT 
decomposition, a dual operator must be constructed. Secondly, a sampled, discretized STFT 
operator partitions an analyzed function's two-dimensional conjugate space into uniform, 
rectangular partitions [19]. In the conjugate space, the spectral bin partition dimensions are 
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inflexibly dependent upon the window function w(t) and to not vary with either translation or 
modulation. Many processes (for example, biological processes) can be characterized by 
components whose bandwidths increase with frequency. Consequently, (4.2) provides 
representation which is less than optimum processes comprised of spectral components of 
varying bandwidths [20], 

For the reasons described above, alternative methods of signal decomposition have been 
suggested. Instead of a representation vector of the form (4.3), employment of a representation 
based on a family of functions [21] 

Va.b(t) = |ar 1/2 V (^) (4.4) 

produces a transformation W: R — » R 2 such that 

V{f(t)} = W(a,b) = la|- 1/2 Jl f(t)vj/’(V)dt. (4.5) 

Similar to the STFT, transformations of the form of (4.5) also comprise a distinct class in the 
Field of group theory. Transformations based on scaling and translation of a common 
representation vector comprise the affine group [18]. 

If, as in the case of digital signal processing, it is desired to restrict the transform to a 
lattice of discrete points, the representation vector becomes 

Vm,„(t) = ao m/2 y(ao m t-nb 0 ). (4.6) 

The representation (4.6) results in a transformation R — » Z 2 , where Z is the set of all 

integers, such that 

W{f(t)} = w(m, n) = ao m/2 [^ f(t)vj/*(a^"t-nbo) dt. (4.7) 
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Additionally, (4.6) constitutes a version of (4.4) sampled at a = a 0 ra for a,, > 1 and at b = 
n bo-a^" 1 for b 0 * 0 . To conform to the conventions of octave-band filtering, is typically 

selected as a 0 =2. The selection of b 0 determines whether {\jr m n } ng z constitutes an incomplete 

set, a complete orthogonal set, or a frame. The representations of greatest interest to 
multiresolution signal processing are chosen such that b 0 = 1, providing unit translations with 
respect to sampled data. 

The time-frequency properties of the representation vector n address some of the 
shortcomings sometimes ascribed to the STFT representation. As the scaling integer m 
increases, the representation \jr m n becomes more and more spread in the time domain, and, 
consequently, more concentrated in the frequency domain. Consequently, projection on a highly 
dilated vector function will provide poor time resolution but sharp frequency resolution. 
Decreasing m causes the reverse effect: Concentration in time and spreading in the frequency 
domain. In the case of a highly contracted representation vector, the transform operator provides 
sharp temporal resolution but poorer spectral resolution. Furthermore, the spreading effect of 
the representation V|/ m n occurs in a logarithmic manner. The bandwidth of the representation 

\y m n be proportional to its center frequency. 
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Figure 4. /—Venn diagram illustration of concept of embedded vector spaces. After [24], 



B. THEORY OF MIJLTIRESOLUTION ANALYSTS 

In order to characterize the vector spaces consisting of the span of {t|/ m n } m neZJ it is first 

necessary to consider another set of basis functions { n } ns z w hich spans another set of vector 
spaces {V m } nsZ e L : (R). The operator A m ., is defined as the projection of some function f(t) on 

V«: 

A m -,{f(t)}=l(V, 4> m . n ) 4>m.n(t). (4.8) 

To be a multiresolution analysis , the set of operators { A m } must satisfy six properties [22,23]: 

1 . A m _, is a linear operator which uniquely and completely approximates f(t) at a resolution 
of 2 m . Consequently, the approximation 

(4.9) 

In words, A m _,{f(t)} contains all of the information about f(t) which can be obtained at a 
resolution 2 m . Repeated projection upon V m does not add or subtract any information to 
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2. Of all possible functions which exist at a resolution of 2 m , A^., (f(t)} is the function 
which most closely resembles f(t): 

v g(t) « V„, ||s(t) - f(t)|| s l|A m .|{f(t)} ■ f(t)||. (4.10) 

3. Approximation of some signal f(t) at one resolution 2 m contains all information 
necessary to approximate it at the next resolution 2 m ~'. This concept suggests a family of 
embedded, closed subspaces: 

v m^. c v m 6 L : ( R ) V m s Z. (4.11) 




Figure 4.2 - Spectral illustration of concept of embedded vector spaces. The vector spaces represented are 
based on Daubechies' orthononnal scaling function on [0, 1 1]. After [25]. 

Illustration of the concept of embedded spaces can be accomplished by either of two 
methods. The first illustration is via Venn diagram. In Figure 4.1, the embedded ellipses 
illustrate two related vector spaces. The outer ellipse represents the span of the vector space V m 
and the inner ellipse the span of vector space V m+1 . As indicated by the diagram, the information 
which can be extracted from projection on vector space V m , t is less than what can be extracted 
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from projection on vector space V m . The two spaces differ by the information lost when some 
function f(t) is approximated at 2 m *' instead of at 2 m . The vector space V m possesses a greater 
density than the vector space V m+1 . Consequently, the approximation at 2 m will contain a greater 
quantity of information about the original function than the approximation at 2 m '" 1 . 

Figure 4.2 provides an illustration of the concept of embedded subspaces from the 
perspective of partitioning of the frequency spectrum. The lower half-band represents the vector 
space V m while the lower fourth-band represents V mH . Approximation A m .,{f(t)} at resolution 
2 m , therefore, entails an approximation based on the spectral components of f(t) contained in the 
lower half of the frequency spectrum below the Nyquist frequency. Similarly, approximation at 
a resolution 2 mH entails a representation based on the lower fourth of the frequency spectrum 
below the Nyquist frequency. Consequently, an approximation of f(t) based on Figure 4.2 at 
resolution 2 m contains only the spectral content of f(t) in the range [0, Jt/2]. An approximation 
of f(t) at resolution 2 m ~‘ contains the only the spectral content of f(t) in the range [0, Jt/4], 

4. The approximation operation is similar at all resolutions. The spaces of approximated 
functions can, therefore, be derived from one another by scaling each approximated 
function: 

V m € Z, f(t)eV m «, f(2-t) e V m .,. (4.12) 

5. The approximation A m .,{f(t)} of a signal can be characterized by 2 m samples per unit 

interval. When f(t) is translated by an amount proportional to 2 ™, { f(t) } is translated 
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by the same amount and is characterized by translations of the vector space projections. 
More simply, 

V n e Z. f(t)<= V m <=> f(t-2‘ m n) € V m . (4.13) 

Equations (4. 12) and (4.13) suggest a family of basis functions <f) mjl (t) similar to that 
characterized by (4.6): 

<t>m.n(t) = 2~ mZ <{)(2 -m t — n) V m, n € Z 2 . (4.14) 

The vector space {<J> m n } ng z e V m consists of integer translates of (The t-dependence has 
been suppressed for compactness of notation.) The approximation operator (fft) } is simply 
a projection of f(t) in the space of the vectors Additionally, (4.12) reinforces the concept 

illustrated in Figure 4.2. Equation (3.7) indicates that a time-domain contraction of a signal 
causes a dilation-or spreading-of that signal's frequency spectrum. The frequency spectrum of 
V m+ „ therefore, occupies half of the bandwidth of the frequency spectrum of V m . Furthermore, 
if the space spanned by V m coincides with a lowpass region of the frequency spectrum, V m _, 
will also coincide with a lowpass region. 

Finally, one additional property remains to complete the definition of a multiresolution 
analysis: 

6. A continuous function f(t) can be initially considered to be represented with infinite 

resolution. Regardless of the scale, all information regarding f(t) is originally assumed to 
be known. Applying the approximation A m { f(t) } results in some loss of information. 
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Increasing the degree of coarseness of A m increases the amount of lost information. 
Consequently, 

lim V m = 0 V k =L 2 (R) (4.15a) 

m— k=~°° 

and 

Jirn V "> = k ^ v k = { 0 }. (4.15b) 

Equations (4.15a) and (4.15b) follow, by induction, directly as consequences of (4.1 1). 
Equation (4.15a) states that, because of (4.1 1), the manner in which {V m } mg z are related to 
each other, lim V m consists of the union of the spans of all { V m } m „ 7 . By induction from 

m— » — » v in ^ m e A J 

(4.11), all vector spaces {v} z are subsets of lim V m . Using the notation of set theory, 
the concept behind (4. 15a) is expressed as 

(V,),,z <= Jta.Vm. 

Consequently, adding to (4.15a) the concept expressed by (4.12), as the resolution of an 

approximation A m ,, the projection on vector space V m , becomes infinite, it is possible to 

represent f(t) with arbitrary precision. Furthermore, employing a concept from Chapter II, 
lim V m is dense in L 2 (R). In other words, given any vector 0. it is possible to find another 

vector <t> m-k * a which is arbitrarily close to 0 mjl . Contained within the union of all definable vector 
spaces {VJ is all the information which is known about a function. 

The indication of equation (4.15b) exactly opposite of that of (4.15a). By induction 
from (4.1 1), lim V m constitutes the least common subset of all of the vector spaces 

m — 

{ V m } me z-C° nse q uen tly, as the coarseness of the approximation A^t)} becomes infinite, all 
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information about f(t) is lost. From the perspective of Figure 4.2, the portion of the frequency 
spectrum spanned by lim V m is contained in the closed interval 0, lim . 

m— *«> |_ m— ► ^ J 

Consequently, at the limit, the approximation of f(t) is characterized by only its DC component, 
or equivalently, lim A m {f(t)} becomes a constant-valued function. 

m — * eo 

Regarding multiresolution analyses, Mallat [22] proved a theorem which provides the 
theoretical foundation for all further development. The theorem states that, given a 
multiresolution approximation of L 2 (R) by projection on vector space V m , there exists a unique 

function <|)(x)€ L 2 (R), called a scaling function, such that {4> mjl (x)} ngZ as defined by (4.14) 
constitutes an orthonormal basis of V . That is, 

m 5 

^$m.n> ^m.k ^ = §k,n ^ n, k € Z. (4.16) 

Furthermore, V m consists of the closed, linear span of n } ne 2 .- 




Figure •/.i—Superimposed plots of Daubechies' orthononnal scaling function <J> 0 0 (t) supported on [0, 1 1] and 

of <t>.,. 0 (t). 
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To characterize the relationship between the scaling function spaces {V m }, it is useful to 
examine the projections of their vectors onto each other. From (4.1 1), (4.12) and (4.13), it is 
evident that 

4> m -l. n ( t ) e v m *. <=* v m (4.17) 

or equivalently, <j> mH n lies within the span of {<J> m „}• Figure 4.3 provides an illustration of a 
scaling function at two adjacent scales. The example presented is based on Daubechies’ 
orthonormal scaling function supported on [0, 1 1]. <+>., 0 (t) is obviously a contracted version of 
4> 0 0 (t) * ts amplitude increased by a factor of Jl . Because of (4. 17), the Fourier series 
theorem (2.17) can be applied to obtain 

4*m+l.n(t) = 2 ^m+1. n > ^m, k ^ k(0 • (4 .18) 

Substituting (4.14), the definition of <t> m n , into the inner product term of the summation (4.18) 
and applying a change of variables to evaluate the resulting integral produces [22] 

^tj)m+i,n, 4*m, k ^ = — ' ^4* 1,0» 4*0. k-2 n ^ • (4.19) 

Equation (4.19) clearly indicates that the coefficients of the series (4.18) are independent of 
scale m. The summation (4.18), therefore, becomes 

4>m+l.n(0 = ^=* 2^4*1. 0> 4*0,k-2 n j ' 4*m.k(t)- (4.20) 

Next, substituting (4.14) for the appropriate terms on each side of (4.20) results in 

2“(m+l)/2 ^(2-tm+l) t - n) = -L. . 2~ ml I (<]>,. 0, 4*0, k-2 n ] ' 4*(2" m t - k) . (4.21) 

Applying the change of variables u=2" m • t-2 ■ n and the translation of indices k=k'-2-n to (4.21) 
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finally yields 



4>(u) = Z ( <t>1.0, 9o.k-2n J<W 2 u-k). 

For later convenience, (4.22) can be rewritten as 

7 <t>(y) = Z h(k)0(u-k) where h(k)=Y^<t>i. 0 , <J> 0 .k j. 



(4.22) 



(4.23) 



Selecting this normalization for h(k) forces Z h(k) = l . Asa result, 



Zh(k) e jt0 

k 



co = 0 



= H(eJ“) 



co = 0 



= l 



Furthermore, since <J)(u) is real, {h(k)} ks z is real. And finally, if <f>(u) is compactly supported, 
then the filter h(k) is a FIR filter. This occurs since, for compactly supported 4>, in the inner 
product term of (4.23), there will be only a finite number of translations which will be evaluated 
as non-zero. The details of the development of (4.19) - (4.23) are presented in Appendix A. 



Equations of the form (4.23) comprise in literature a class of equations referred to as 
two-scale difference equations [26] or as dilation equations [21]. This class of equations will be 
considered in greater detail in Chapter V within the context of basis functions for multiresolution 
analyses. From the preceding development it may be concluded that the functions <|) which 
satisfy the properties for a multiresolution analysis consist of solutions to two-scale difference 
equations. 

To obtain a multiresolution transform it is necessary to combine (4.20) with (4.23) to 

obtain 
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(4.24) 



<JWn(0= v 2 Ih(k-2n)<J> m . k (t). 

k 

Applying Parseval's equality (2.18) to (4.24) results in a relationship between the projections of 
some function f(t) onto V m and V m ., [22]: 

(f, Om+i.n j = J2 I h(k-2n) (f, 0 m . k j. (4.25) 

By (4.25), given A, n .,{f}, the determination of the coefficient (f, <5 m . t n ) for each term of the 



series for A^Jf}, requires filtering and decimation of the sequence { (f, (f> m n )} ngZ . The block 

diagram for the operation (4.25) is depicted in Figure 4.4. The similarity between Figure 4.4 
and a channel of a QMF analysis bank is evident. 



{U4>m,n)}„ 




)}. 



m+ 1 ’ 11 rt€ Z 



Figure 4.V--Block diagram of multiresolution transfonn characterized by (4.25). 

Now, to complete the development of the concept of multiresolution analysis it is first 

necessary to consider some of the consequences of the preceding development. First, consider 
the concept of embedded subspaces characterized by (4.1 1). Let W m be the orthogonal 



complement of V m such that 



v. nw, = 0 


(4.26a) 


and 




V ®W = V ,. 

m m in- 1 


(4.26b) 



From (2.26b), by induction, any scaling function vector space V, can be decomposed as [23] 

Vj = ® W m . (4.27) 

m= J+ 1 
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And finally, the direct sum of all definable {W m } m€Z is dense in L 2 (R) [21]: 



® W m =L 2 (R). (4.28) 

m=-«» 

Figure 4.5 presents a two-scale example of (4.26b) from the perspective of spectral 
contents of V m and W m . The complementary subspace W m is defined such that it is closed and is 
spanned by a set of vectors {V m , n }nez- Members of the set of vectors (Vm.nl are related to each 
other in a manner similar to (4.14): 

Vm,n(t) =2" nv2 \|/(2" m t-n) V m, neZ 2 (4.29) 

Overlap exists between the spectral regions characterizing V m+ , and W m+1 . This overlap occurs 
since the plots consist of frequency plots of FIR filters. However the FIR filters represented are 
related to orthogonal basis function and consequently, the filters for V m+1 and W m+1 are 
themselves orthogonal to each other in the time domain. 




Digital frequency, (multiple of f ) 

Figure ¥.5~Spectral content of example of concept of orthogonally complementary, embedded subspaces. The 
example is based upon Daubechies' wavelet and scaling function on [0, 1 1]. The lowpass processes are 

bn, n and <tWn and the bandpass process is tiv, „ 
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The set of vectors {Vm.JneZ constitutes the set of wavelet functions for the 
multiresolution analysis. From the relationships between the spaces W m of wavelet vectors n 
and the spaces V m of scaling function vectors 4> m n a number of relationships between the sets of 

vectors <t> m n and \jr m n may be deduced. Between the scaling and wavelet function vectors, 
follows from (4.26a), the relationship exists: 

(4V..V.J-0. (4.30) 

Within a common scale, therefore, wavelets and scaling functions are orthogonal to each other. 
Furthermore, as a consequence of (4.1 1) and (4.26b), 

0„,„, ¥,.*) (4.31) 

In words, unlike scaling functions which are only orthogonal to each other with respect to 
translation within a common vector space V m , wavelets are orthogonal to each other with respect 
to translation and scale. Wavelets not contained within the same vector space W m are by 
definition orthogonal to each other. Finally, because \|/ m n € W m c V m . p 

(<t> ra ,n> V k < m - 1. 

More generally, if the wavelet V|/ k ■ is supported on the interval [0, L] e R, then, 

(4> ra , n . vj * 0. V k < m- 1 , I n - j | < L. (4.32) 

Wavelet vectors lie within the span of higher resolution spaces of scaling functions. 
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The projection of some signal f(t) on the wavelet space W m constitutes the detail signal 
[22] and contains the difference of the information contained in A m .,{f} and that in A m {f}. The 
detail signal at resolution m is denoted by 

D m -l{f(t)}=I ( 7 , M/m.n jvm.n(t). (4.33) 

Equation (4.26b) implies that [21] 



D m {f} = A m . I {f}-A m {f}, 



(4.34a) 



or equivalently, 



I 




Vm+l.n 




(4.34b) 



Substituting (4.24) into (4.34b) and expanding the result produces a time-domain form for the 
wavelet decomposition: 






4>m.n, f )<]> m , n -2 Ih(k-2n) h(j — 2 n) ( <p m . k , f J 4> 



m.j 



. (4.35) 



Obviously, except for simple cases such as the Haar basis (demonstrated in [21]), (4.35) is not 
easily separable. In order to obtain a relationship between scaling and wavelet functions it 
therefore becomes necessary to examine their Fourier-domain properties. 



Equation (4.22) represents the starting point for considering the frequency-domain 
properties of the scaling function <{). If the Fourier transform of (f> is defined as in (4.1) 

0(W)=!1 0(t)e" joK dt, 



evaluating the Fourier transform of each side of (4.23) produces [23] 

4>(2co) = H(e i<0 )<J>(co) (4.36b) 
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where the Fourier series H(e iu) ) of h(k) is defined as 

H(e j “)=Ih(k)e-J° k . (4.37) 

k 

Next, Poisson's summation can be employed to present an alternative expression for the 
orthogonality of 0 with respect to translation [23], This development begins with the complex 
Fourier series expansion for a train of Dirac impulse functions 5(co) [28]: 

X5(co-27ik) = 7 -Se jkw . (4.38) 

k 2K k 



Poisson's summation involves convolving a function with each side of (4.38). Convolving the 



expression 



<t>(co) 



with the right-hand side of (4.38) produces 



<J>(co) 



*E5(a>-2jck) = I 

k k 



<t>(co-2 Jtk) 



(4.39) 



Furthermore, convolving 



<Kco) 



with the right-hand side of (4.38) results in the expression 



i>(co) 2 * £ eJ k “ = 1/ 0(5) 2 e ik( ^> d^ 



= £ej k(0 ||| 0(u)0(v)e j(v u) ^ dudv e j k ^ d^ 

k 



(4.40) 



Applying the change of variables of integration w = v - u and regrouping the terms contained in 
the second line of (4.40) produces 

<t>(co) * £ej kft) =£ e jkw }J 0(u)0(u + w) J ej (w_k) ^ d^ dudw. (4.41) 

k k 

Now, since 

| e j(w_k) ^ d^ = 27t5(w-k), 

by the integral sifting property of the Dirac delta function 5(w), (4.41) becomes 



0(co) *Ie jk<d = 2^1 ej kw If 0(u)<))(u + w) du 5(w-k) dw 



= 2;rX ej k(0 | 0(u) 0(u + k) du 

k 



(4.42) 
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Finally, by (4.16), the integral in the final line of (4.42) is simply (<{> 0 0 , <j> 0 J = 8 0 k where 5 m n is 
Kronecker's delta function. Consequently, the series in the final line of (4.42) contains only one 
non-zero term: when the index k = 0. Therefore, 



<t>(w) 



*IeJ ku = 27t. 

k 



(4.43) 



Substitution of (4.39) and (4.43) back into (4.38) produces 

2 



x 

k 



<]>(CD-2 7tk) 

Applying the results of (4.44) to (4.36) yields 



(4.44) 



I 

k 



4>(2co-2jrk) = X Hle^' 110 ) 4>(co-7tk) 

k 1 1 



= l . 



(4.45) 



The right-hand side of (4.45) can be expressed as 



X 

k 



H(e J 



j (g>- 2 k rc) 



) 



<j>(co - 2 ktr) 



+ X 

k 



H(e J 



j (co-{2k+l)rc) 



) 



(j)(CL) -(2k+ 1 )ji) 



= 1. 
(4.46) 



Since, as a consequence of its definition, H(e i “) is periodic with respect to 2-n, the terms of 
| p^uo-’Kk)) 1 2 an( j | ^ e j(o>-( 2 k+nit^ | : can f actorec i out 0 f t h e se ri es on the left-hand side of 

(4.46) producing 

I H( e j«) I 2 X I <j>(co- 2 k 7 t) 2 +| H(e iM ) 2 X $(co -(2k+ 1 )jc) 2 = 1. (4.47) 

1 1 k 1 1 k 

The remaining two sums in (4.47) are identical to (4.44) and are evaluated as unity. Therefore, 

(4.47) becomes 

I H(ei“) 2 + | H(e Jl( °-^) | 2 = 1. (4.48) 
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Comparing (4.48) with (3.20b) reveals that the frequency response H(e ,<tf ) for the 
approximation filter defined by (4.23) is power-complementary with a version of itself shifted in 



the frequency domain by K. Therefore, Hfe 1 C) ) either is constant or possesses a half-band 
frequency response. Furthermore, if the scaling function (p is normalized such that 

7 



m 



1 , then, as a consequence of (4.36b), 

I H(e J ") | 2 = 1. 



(4.49) 



C. MALLAT ’S ALGORITHM FOR MULTIRESOLUTION ANALYSIS 

It remains to characterize the wavelet functions. This can be accomplished by 

completing the development of the multiresolution decomposition begun in the previous section. 
The first true multiresolution algorithm was developed by Mallat [22], Two other techniques 
will be considered in the next section which were not directly derived from the mathematical 
definition of a multiresolution analysis (4.8)-(4. 13). 

The multiresolution decomposition process begins with the projection of a function f(t) 
on the scaling function space of finest resolution {4> 0 ,„}„ 6 z [21]. The functional expansion 
appears as 

f(t) = A_i{f(t)} = Sco.n<>(t-n) where c 0 . n =^f, j. (4.50) 

Generally, the decomposition is performed on a sampled data sequence and the scale of the 
initial approximation is defined such that the first equality in (4.50) holds. Since the basis 
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vectors (|) 0 n € V 0 , the original approximation (4.50) is also clearly contained within the span of 

v 0 - 

The next step in the multiresolution decomposition requires the division of the scaling function 
space V 0 into orthogonally complementary subspaces V, © W, = V 0 . Re-arranging (4.34a) leads 
to 

A,{f(t)}=A 0 {f(t)} + D 0 {f(t)}. 

Consistent with the definition (4.8) of the approximation transform A m and the definition (4.33) 
for the detail transform D m , 

Ao{f} - SCi.k^i.ic 

k 

and 

Di){f} = Zbi > i t \|/'i.i t 

k 

where c, k = (f, <J>j J and b, k = (f, % j. A relationship between q>, k and (f> 0 n was presented in 
(4.23). Furthermore, (4.32) showed that V|/ 1>k lies within the span of V 0 . Consequently, 
following development analogous to that for (4.23), the Fourier series expression for wavelet 
function becomes, 

y\jr(y) = Ig(k)<t>(u-k) where g(k)=y (yi.o, <Kk), (4.51a) 

since 

^M^rrv+l,n> $m. k ^ = — 1.0 > $0, k-2o ^ • (4.5 lb) 

Since, if \y m n is real, the filter coefficients {g(k)} kg z are strictly real. Additionally, Fourier 
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transformation of (4.51a) produces a relationship analogous to (4.36): 

y(2co) = G(e J “)<i>(co). (4.51c) 

Finally, the form of (4.24) for the wavelet function appears as 

V«tti.»(t) = J2 Ig(k-2n)<t> m . k (t). (4.52) 

k 

Applying Parseval's equality (2.18) to (4.24) and (4.52) in order to calculate c, Q and n in terms 
of c 0 „ produces 

C].„ = yj • Eh(k-2n)co k 

k 

and 

b i,n = 72 •,Ig(k-2n)c 0 .k . 

k 

The approximation based on the expansion of the terms {c, „•<)>, n } constitutes a "blurred" version 
of the expansion of the terms (c 0 n -(f> 0 n }. The detail which is lost in the approximation A 0 { f} is 
contained in an expansion of the terms {b, n }. This process can be repeated as many times 

as desired. The resultant general expression becomes 

Cm+i.n = 72 •Ih(k-2n)c m . k 

k 

and , (4.53) 

bm+l,n = ■ X g(k — 2 n) C m , k 

k 

where each successive iteration generates in an approximation of f(t) containing half the 
resolution of the previous approximation. At each iteration, the approximation and detail of the 
signal are related to the previous approximation by 

A^,{f}-AJf} + D„,{f}. (4.54) 
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Furthermore, the operations described by (4.53) closely resemble the analysis bank for the QMF 
bank structure investigated in the previous chapter. The approximation and detail operations at 
each resolution consist of filtering a common sequence with two different filters. The output of 
each filter is then decimated. The resemblance of QMF structure to the operation characterized 
by (4.53) cannot yet be asserted to be exact since the appropriate relationship between the filters 
whose impulse responses are h(k) and g(k) yet to be established. 

To reconstruct a signal from its decomposition, the definition (4.8) of the approximation 
operator is employed. The expansion coefficient c m k is, by definition, 

c„. k - (A m ,{f), O = (A„{f>, + (DJf), (4.55) 

Substituting into (4.55) the expressions for the series A,,, and D m produces 

Cm.n = £ Cm+l.k ^m,ni ^m+l.k ^ + Z b m+ |, k n> Vm+1. k ^ • (4.56) 

The inner product terms in (4.56) are precisely those encountered in (4.19) and in (4.51b) except 
that the positions of the translation indices n and k have been reversed. Consequently, (4.56) is 
equivalent to 

c m ,n = ^=r-£h(n-2k)c m+ltk +-=-• Sg(n-2k)b m+ i. k (4.57) 

J2 k ,’2 k 

The operation in (4.57) entails expansion of the detail and approximation series coefficients 
followed by the application of filters whose impulse responses are the time reversals of those 
used in (4.53) for signal decomposition. Again, a resemblance to the QMF bank synthesis 
structure is noteworthy. 

Examining the decomposition operation in the frequency domain provides further insight 
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into the nature of the wavelet functions [23], First the Fourier transform A 0 (w) of the 
approximation 



A l (t)=A. 1 { fift) }=f(t) is expressed by 

A-i(CO) = Zco.n jr. <t>0.n(t)e -J<ot dt 

= Z Co n e" J<on f“ 6(t)e - j wt dt (4.58) 

n ' 

= C 0 (eJ a )$(co) 

A A A 

where C m (e JC0 )=Sc m k e _,tok . Applying similar analysis, A 0 (co) and D 0 (co), the Fourier 

k 

transforms of the approximation A 0 (t)=A 0 {f(t)} and D 0 (t)=D 0 {f(t)}, respectively, appear as 



A 0 (co) = Zd.k <t>i. k (t)e j at dt 

k 

= -4rZc 1<k j:^ <KT-k)e-*“‘dt 

k 

= J2 Zci,i c e“» 2 “ k £. 4>(u) e"i“ 2u du 

k 

Ao(co) = yi C,(ej 2 “)i(2co) 

and 

Do(co) = J2 B,(eJ 2ffl )v(2co) • 
where B m (e J “)=Zb m k e“ jak . By (4.54) therefore, it follows that, 

k 

A A A 

A_i (co) = A 0 (co) + D 0 (co) , 

or equivalently 



f(co) = C 0 (eJ “) J(co) = /2 C ! (e^ 2 “) ^(2 to) + 71 B , (eJ 2 “) v(2 w) . 
Substitution of (4.36) and (4.51c) into (4.60) produces 

C 0 (eJ“)<t>(G)) = J2 C,(eJ 2 “)H(ej“)^(co)+ VI B,(ei 2u )G(eJ a )<j>(co), 

A 

or, after dropping the common factor of 4>(oo), 

C 0 (ei °)-Jl C,(eJ 2 “)H(eJ (0 )+ Jl B,(ej 2(0 ) G(ei“). 



(4.59a) 



(4.59b) 



(4.60) 



(4.61a) 
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From their definitions, C 0 (e ,a> ) e L 2 ([0, 2-rc]) and is periodic with respect to 2 rc while 
C/e 1 2 “), B,(e )2 “) € L 2 ([0, 7t]) and are periodic with respect to n. Therefore, it follows from 
(4.61a) that 



Co(e J<0 ) 


= J2 


H(ei“) 


G(eJ“) 


’ Ci(e j2<0 ) " 


_ C 0 (e- i(< ‘ >+ ' u ) 




H(e J(6Wl >) 


G(eJ (0 *-' t) ) 


B,(e j2<i> ) 



From (4.61b), define 



(4.61b) 



H(e J “)= 



H(eJ“) Gfe^) 
H( e i (0>+ '' ) ) G(e J(0>+rt) ) 



Now, if the Hie*") is unitary in the sense described in Chapter III, then multiplication of each 



side of (4.61) with its Hermitian transpose yields 



C 0 (e J(D ) | 2 + 



C 0 (eJ (a>Ht) ) I 2 = 2 




(4.62) 



Integrating (4.62) over the interval [0, rr] produces 

jo | Co(eJ“) | 2 +| C 0 (eJ (owt *) | 2 dco = 2 | C,(ej 2 “) | 2 +| Bi(e j2 “) | 2 dco 

which is equivalent to 

Jo" | C 0 (e Jl °) | 2 dco = 2 Jg | C,(ej 2 “) | 2 dco-i-2 | B,(eJ 2<0 ) | 2 dco. (4.63) 
Equation (4.63) is an expression of the orthogonality of the decomposition: The energy in the 
approximation Ao(f(t)} is simply equal to the sum of the energies in the approximation A,{f(t)} 
and the detail D,{f(t)}. No energy is present in terms comprised of the cross-product 

I C,(e )2 “) | • | B 1 (e r2 “) | . Furthermore, the specification of H^”) as unitary is exactly 



equivalent to the specification that the alias compensation matrix for the two-channel QMF bank 
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is exactly equivalent to the specification that the alias compensation matrix for the two-channel 
QMF bank be paraunitary. Consequently, a single-stage decomposition followed by a 
single-stage reconstruction employing Mallat's algorithm is exactly equivalent to a QMF bank. 

Moreover, the paraunitary nature of H(e' w ) implies 

| | ! + | G(e'“> | ; = 1 (464) 

H(e J “)G(e _J “) + H(e ,(t>_, ' , )G(e' J(0_ ’ t ') = 0 

One solution for Gie 1 a ) as constrained by (4.64) is given by 

G(e ico ) = e‘ JCdL H(e )to_,t) ) (4.65) 

where L represents the length of the FIR filter h(n). Selection of (4.65) dictates that the filter 
whose impulse response is g(n) will be a FIR filter of length L and that, by (4.51c), the wavelet 
function \j/(t) will have exactly the region of support of the scaling function cp(t). 

D. EARLY MULTIRESOLUTION ALGORITHMS 




Figure 4.6 — Block diagram of process by which detail is extracted from signal through the use of the Laplacian 

pyramid algorithm. After [29], 

In this section, two early multiresolution techniques— the Laplacian pyramid and the a 
trous algorithm— will be considered. The first, the Laplacian pyramid, represents an algorithm 
developed in the early 1980’s for the purposes of image coding. Conceptually, the Laplican 
pyramid very closely resembles Mallat's multiresolution algorithm presented in the previous 
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section. In fact, Daubechies credits the Laplacian pyramid with providing Mallat with 
significant for the development of his own multiresolution algorithm [21]. 

The Laplacian pyramid, as with Mallat's multiresolution algorithm, entails iterative 
filtering of a sequence to progressively smooth out the rapidly varying components [29], At any 
stage, the next coarser approximation is obtained from 

fit + i(n) = Iw(2n-m)fk(m). (4.66) 

m 

Unlike Mallat's algorithm, however, where the detail is extracted through filtering, the Laplacian 
pyramid extracts the detail through re-expanding and filtering the approximation and subtracting 
the result from the previous approximation. This process is illustrated in Figure 4.6. 

Mathematically, extraction of the detail lost A k+1 (n) when f k (n) is approximated as f k+1 (n) is 
expressed as 

A k+l (n) = f k (n)-2I w(n -2m) f k+l (m). (4.67) 

m 

The series expression on the right-hand side of (4.67) is recognizable as the same form of 
equation as appeared in (4.57). The form represents expansion of a sequence followed by FIR 
filtering. The factor of two is introduced to cancel the factor of 1/2 introduced by the 
decimation operation as indicated by (3.6). As with Mallat's algorithm, iteratively repeating, 
results in the cascaded tree structure of Figure 4.7a. Reconstruction of a decomposed signal is 
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Figure 4. 7a~Block diagram illustrating three-stage decomposition by Laplacian pyramid. 




Figure 4.7b— Block diagram of structure for reconstructing sequence decomposed by three-level decomposition 

depicted by block diagram of Figure 4.7. 




Figure 4.£~Lattice illustrating concept behind generation of kernel w(n) used for Laplacian pyramid 

decomposition and reconstruction. After [30]. 
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accomplished by simply reversing (4.67) and reintroducing to f kM (n) the detail A kH (n) removed 
during decomposition: 

f k (n) = A k+ i(n) + 2- Sw(n-2- m)f k+ i(m). (4.68) 

m 

Figure 4.7b illustrates the process by which a sequence decomposed by Figure 4.7a is 
reconstructed. 

Figure 4.8 provides a graphical illustration of the averaging process employed in the 
Laplacian pyramid. The Laplacian pyramid decomposition involves the approximation of some 

sequence f 0 (n) by an averaged version f,(n). The technique developed by Burt and Adelson 
[30], in calculating the values of the nodes in f,(n), endeavors to consider each node in f 0 (n) with 
an equal weight. In the case of a five-point weighting sequence w(n), each node in fj(n) is 
calculated from an average of its five nearest neighbors in f 0 (n). For example, the value of node 
f 0 (n) in Figure 4. 1 8 is evaluated as 

fj(0) = w(-2)-f 0 (-2) + w(-l)-f 0 (-l) + w(0)-f o (0) + w(l)-f 0 (l) + w(2)-f 0 (2). 

Additionally, obtaining f,(-l) and f,( 1 ), requires evaluation of 

fj(-l) = w(-2)-f 0 (-4) + w(-l)-f 0 (-3) + w(0)-f 0 (-2) + w(l)f 0 (-l) + w(2)-f o (0) 

and 

f,(l) = w(-2)f 0 (0) + w(-l)f 0 (l) +w(0)f 0 (2) +w(l)f 0 (3) + w(2)f 0 (4). 

Now f o (0), an even-numbered node, contributes to the computation of three nodes in f,(n): to 
f,(0), to which it is directly adjacent across scale, and to f , ( - 1 ) and f,(l). However, f 0 (l), an 
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odd-numbered node, only contributes to the contribution of two-nodes in f,(n): to f,(0) and to 

f,d). 

As stated previously, the objective is to define {w(n)} n such that the total weighted 
contribution by each node in f 0 (n) to the computation of nodes in f,(n) is equal. Therefore, the 
requirement that the total weighting factor for all contributions by f 0 (0) to equals the total 
weighting factor for all contributions by f 0 (l) implies that 

w(- 1 ) + w( 1 ) = w(0) + 2-w(2). (4.69a) 

The left-hand side of (4.69a) represents the sum of the weighting factors for all contributions by 

f 0 (l) to all nodes to which it contributes in f,(n). Similarly, the right-hand side of (4.69a) 
indicates the sum of the weighting factors for all contributions by f 0 (0) to all nodes to which it 
contributes in f,(n). Inductively, for a five-point weighting sequence the sum of the weights for 
all contributions by any even-numbered node in f 0 (n) will equal the right-hand side of (4.69a) 

while the sum of the weights for all contributions by any odd-numbered node in f 0 (n) will equal 
the left-hand side of (4.69a). 

An additional constraint was imposed on the selection of the weighting sequence w(n): 

I w(k) = 1 . (4.69b) 

k = -2 

The constraint imposed by (4.69b) ensures that, for any node in f,(n), the sum of the weights by 
all contributing nodes in f 0 (n) will be unity. Consequently, the average energy in f 0 (n) is 
preserved by the approximation f,(n). 
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(4.70) 



To satisfy (4.69a) and (4.69b), Burt and Adelson selected 




Each odd-numbered node contributes twice, each time with a weighting of 1/4. The 
even-numbered nodes contribute with a weighting of "a" when directly adjacent across scale 
while contributing twice with a weighting of 1/4 - a/2 when separated by scale and translation. 
Consequently, each node contributes with a total weight of 1/2. 

Considering the top branch of the structure in Figure 4.6 from a two-sided Z-transform 
perspective provides insight into the frequency-domain character of the Laplacian pyramid. The 
transform of the weighting sequence can be expressed as 

W(z) = ^-yj(z 2 + z- 2 ) + |(z + z- 1 ) + a. (4.71) 

Applying (3.6), decimating (4.71) produces 

W i2 (z) = y [W(z 1/2 ) + W(-z I/2 )] 

= )( z + z ~') + a 

which after expansion becomes 

[W i2 ] T2 (z) = (j--f ) (z 2 +z" 2 )+ a. (4.72) 

The effective transfer function T(z) of the top branch of the structure in Figure 4.6 can be 
expressed as 

T(z) = 2-[W i: ] T ,(z)-W(z). (4.73) 

Through trial and error, Burt and Adelson observed that for the choice of a = 3/10, repeated 
applications of the transfer function (4.73) produced an impulse response with approximately the 
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same shape with each iteration and which ultimately, approached a Gaussian shape. For this 
choice of "a", the system transfer function becomes 

JL Z 8 + X Z 7 + -L z 6 + l Z 5 + II Z ' + 1 Z 3 + 1 L z : + _L Z+ JL 

$0 20 25 5 SO S 25 20 SO 

Furthermore, from Figure 4.6, by inspection, the expression for the detail sequence A kM (n) is, in 
terms of f k (n) and t(n), 



A k+I (n) = 1 t(n-k)f k (k)-f k (n). 

k 

Equivalently, if the transfer function for the entire structure of Figure 4.6 is defined such that its 
impulse response is d(n), then d(n) = t(n) - 5(n) where 5(n) is Kronecker's delta function. 
Consequently, the transfer function D(z) becomes 



X z s + X z 7 + ^. z 6 + i z 5_Ji z 4 + i z 3 + x z : + _L z+ JL 
__ so :o 25 s so s :s :o so 

From this analysis follows a conceptually clearer but computationally less efficient structural 
equivalent for Figure 4.6. This equivalent structure is indicated in Figure 4.9. 




Figure 4.9-Equivalent structure to that shown in Figure 4.6. 

As a method of time-frequency decomposition, the Laplacian pyramid provides 

performance for processes containing high-frequency components which inferior to that which 
shall later be observed for Mallat's algorithm. Figure 4.10 illustrates the partitioning of the 
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frequency spectrum which occurs when a signal is decomposed by the Laplacian pyramid. The 
spectral bins which are formed do not follow the pattern of an even division of the frequency 




Figure ‘/./£)--Panioning of frequency spectrum perfonned by three-stage Laplacian pyramid of Figure 4.7. 
spectrum at each stage as would be anticipated by operations involving factor-of-two decimation 

and expansion. In fact, approximately three fifths of the frequency spectrum below the Nyquist 
frequency is contained within the highest frequency-bin. Consequently, for classification 
applications, a Laplacian pyramid-based analyzer would provide poor localization of frequencies 
above approximately 0.2 f $ . However, for an over-sampled sequence comprised primarily of 
lower frequency components, a Laplacian pyramid detector would likely provide slightly better 
spectral resolution than would be afforded by a multiresolution detector which divides the 
spectrum evenly at each stage. For a three-stage, even-division multiresolution scheme, the 
region in the frequency domain from 0 to 0.125-f s would be divided into two spectral bins. For 
the scheme of Figure 4. 10, the region in the frequency domain from 0 to approximately 0. 1 5-f s is 
divided into three spectral bins. 
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Figure 4. /./--Block diagram for full expansion by channel approach for reconstruction of a sequence decomposed 
by the Laplacian pyramid. Structure is equivalent to that of Figure 4.7b. 

In general, multiresolution decomposition schemes which involve factor-of-two 
decimation generate lattices which resemble Figure 4.8. At each stage, the approximation 
density is reduced by a factor of two. For information transmission systems, this represent the 
primary advantage of such systems. However, for display of time-scale decompositions 
involving mesh or contour plots, most computer graphical routines require the insertion of zeros 
in order to create a lattice of points of constant density. For mesh plots, such as those provided 
by Matlab, a lattice of the form of Figure 4.8 will appear as rows of fin-like structures of 
varying density. Because of the rapid transition to zero at each lattice nodes, contour plots 
essentially provide binary indications of a sequence's time-scale content: Indication of the 
presence of energy at a node is given, without indication of the relative magnitude of the 
sequence at that node. 



Figure 4.1 1 provides a block diagram of a full expansion by channel scheme for 
preparing the lattice similar to that of Figure 4.8 for plotting. The structure of Figure 4. 1 1 
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represents the result of separating each branch of Figure 4.7b at its summation points. The 
consequence of implementing a structure like that in Figure 4. 1 1 is that sparse rows in a lattice 
like that in Figure 4.8 will be interpolated resulting in a uniform density of lattice nodes. 




Figure 4. 12-Time plot of 256-point test sequence generated by (4.74a) for demonstration of Laplacian Pyramid. 




Figure 4.13 - Plot of power spectral density of 256-point test sequence generated by (4.74a). 
Furthermore, the reconstructed sequence from the structure in Figure 4.7b is 

accomplished by 

fo(n) = If m »(n). 
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A 

In other words, the reconstructed sequence f 0 (n) consists of the sum of the content of all the 

channels f m) (n) obtained from Figure 4.11. 

To demonstrate the operation of the Laplacian pyramid, two test sequences were 
employed. The first sequence was constructed from low-frequency components: 

SLF(n) = ^ 1 - e _n/6J j • e 2n/75 • [2 cos(2 rt n) + sin(2 • n ■ n) + y cos(2 ■ n 4^- n)] . 

( 4.74a) 

A time plot of the sequence is presented in Figure 4.12 and a power spectral density plot appears 
in Figure 4.13. The test sequence was chosen such that the spectral content would lie 
predominantly in the lower portions of the frequency spectrum. Test sequence (4.74a) contained 
one spectral component, which completes only one complete oscillation during the duration of 

the sequence, at tr/128 (= 0.004-f s ). Additional components are present at 9-tr/256 (= 0.01 8-f s ) 
and at 3 1 -7C/256 (= 0.06-f s ). 

The second sequence, identical to (3.47) employed in Chapter III, Section C is repeated 
here for convenience: 

SHF(n) = ^1 -e _n/64 j • e 2n/75 • [cos(2 • n ■ • n) + j cos(2 • n ■ ■ n) + j cos(2 • k ■ ■ n)] 

(4.74b) 

A time plot of (4.74b) is presented in Figure 3.23 and its power spectral density is plotted in 
Figure 3.24. Test sequence (4.74b) contains spectral components at 2-7T-27/256 (= 0.105-f s ), 
2-K-55-K/256 (= 0.215-Q and at 2-tr-l 13-K/256 (= 0.44-f s ). 
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Figure 4.14a presents a surface plot depicting the eight-scale, full expansion by channel 
of the Laplacian pyramid decomposition of the test sequence generated by (4.74a). Despite the 
fact that the spectral components of the lowpass test sequence (4.74a) are quite close together, 
the Nevertheless, the Laplacian pyramid is able to distinguish between the three separate 
components. One peak is evident at around sample 32 at a scale of five to six, a series of 
approximately four peaks are evident at a scale of three to four from samples 10 to 60 and a 
series of peaks corresponding to the high-frequency components are apparent at scales one to 
two from samples 20 to 120. 




Qj 

Qj 

-J 



0; 

C 

§ 

a 



Figure 4.74a~Surface plot of full expansion by channel of Laplacian pyramid decomposition of 256-point, 

lowpass test sequence generated by (4.74a). 

Figure 4.14b illustrates the Laplacian pyramid decomposition of the highpass 256-point 
test sequence (4.74b). As discussed within the context of the frequency partition diagram, 



90 



Figure 4.10, the Laplacian pyramid produces significantly poorer resolution for sequences 
comprised of high-frequency components. In fact, for the highpass test sequence (4.74b), the 
two spectral components with the highest frequencies fall within a single spectral bin. 
Consequently, as Figure 4. 14b indicates, it is very difficult using the Laplacian pyramid to 
resolve spectral components residing in the upper two thirds of the frequency spectrum below 
the Nyquist frequency. 
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Qj 
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Figure 4. 14 b ~ Surface plot of full expansion by channel of Laplacian pyramid decomposition of 256-point, 

highpass test sequence generated by (4.74b). 

The reconstruction error for the Laplacian pyramid operation on the test sequence was 
consistently outstanding. The mean square reconstruction error as defined by (3.32) is plotted in 
Figure 4.15. As indicated by Figure 4.15, the reconstruction error never rose above -317 dB. 
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This accuracy roughly corresponds to the numerical precision of Matlab which was used for this 
demonstration. 




Number of Oecomoosition Staaes 



Figure ^.75-Trend of mean square reconstruction error for Laplacian pyramid operating on 256-point lowpass 

test sequence of (4.74a). 

The a trous algorithm represents the first multiresolution analysis technique explicitly 
based on an affine-type representation vector. This technique consists of evaluating a 
discretized approximation of the continuous wavelet transform. Specifically, the starting point 
for the a trous method is the discrete wavelet series [31]: 

w(m, n) = 2 -m/2 Zy*(2 -m k-n)s(k). (4.75) 

k 

Essentially, the discrete wavelet series consists of the discrete correlation between some 
sequence s(n) and translations of a sampled wavelet function 2' m/2 y(2' m k - n) at scale m. 
Evaluating (4.75) for w(l, n) produces 

w(l,n) = -L-S\|/*(4-n)s(k). (4.76) 

In (4.76), it is obvious that \j/(y - n) = y(yp-) . Consequently, (4.76) becomes 

w(l,n) = 4-Iy*( i f ! -)s(k). 

V2 k 2 
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Now, if \j/(k) is known for integer-valued arguments, its values at half-integer arguments 
can be approximated through interpolation. If f(k) represents the impulse response for an 
interpolation filter, then \y(k) at half-integers is approximately 

V(t)» & Zf(2k-n)v(k). (4.77) 

“ k 

where the conjugation and time-reversal of the filter and the factor of Jl have been introduced 
for later convenience. The form of the summation term on the right-hand side of (4.77) is 
recognizable as expansion of a sequence vj/( k) followed by convolution with a filter whose 
impulse response is f(-k). 

To obtain the best approximation, special conditions are imposed on the filter represented 

by f(k). Shensa [31] applies the term a trous filter and employs the notation 

f(2 • k) = -=r5 0 ,k, 

*2 

where 5 0 k is Kronecker’s delta function. Perhaps a slightly more illuminating manner in which 

to express this concept is to use the notation 

fi 2 (k) = -Mo. k . (4.78) 

V 2 

In other words, the decimated impulse response f(k) is non-zero only for k = 0. Equivalently, 
the only non-zero, even-numbered coeffecient for the filter f(k) corresponds to k = 0. The 
condition (4.78) is necessary in order to ensure that known values for the interpolated sequence 
v|/(k) are unaffected by the interpolation process. When evaluating f (2 k - n) \y(k), if n is an 
even integer, an argument for which \j/(n/2) is known, 2 k - n will also be an even integer, and 
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the series (4.77) will contain only one non-zero term, for n = 2k. Consequently, (4.77) will be 
evaluated as 

v(f ) = J2 n o)v(f). 

Hence, the factor of 1/ /2 in (4.78). On the other hand, if n is an odd integer, an argument 

value for which the interpolation must be calculated, the series (4.77) will be evaluated as a 
weighted sum of the function values for \j/( n/2) for surrounding even values of n. The a trous 
condition (4.78) assumes a position of importance in the theory of FIR filters. The condition 
(4.78) is equivalent to the half-band condition. A half-band filter is a symmetric, FIR digital 
filter whose impulse response satisfies (4.78). Mitzner [32] showed that filters satisfying (4.78), 
of necessity, also satisfy 

F(e J “) + F(e^ l ‘ 0_ ' t) ) = 2 f(0) where F(ei u ) = If(k)eJ“ k . (4.79) 

k 

Consequently, the interpolation filter whose impulse response is f(k) is also a half-band filter. 
Furthermore, Shensa showed that the filters for Daubechies' orthonormal scaling functions 
belong to a special class of filters satisfying (4.78) known as Lagrangian interpolation filters. 
Specifically, if c(n) is the impulse response for the filter producing one of Daubechies' 
orthonormal scaling functions, the corresponding Lagrangian interpolation filter is obtained from 

f(n) = -=- • Ic(k)c(k-n). 

,2 k 



Continuing the development of the a trous algorithm, substitution of (4.77) into (4.76) 
produces 



w( 1 , n) = IIf(2j + 2n-k)t|/*(j)s(k). 

J k 



(4.80) 
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Now, if a filter g(n) is defined such that 

g(n) = vj/'(-n), (4.81) 

and a change of indices m = j + n is applied to (4.80), it follows that 

w(l, n) = Ig(n-m) If(2m-k)s(k). (4.82) 

m k 

Equation (4.82) simply describes decimation of s(k) followed by filtering with a FIR filter 
whose impulse response is f(k). The result is then further filtered by an additional FIR filter 
whose impulse response is g(n). Furthermore, the second FIR filter is the conjugated 
time-reversal of the sampled wavelet function. Therefore, filtering with g(n) is equivalent to 
evaluating the correlation with y(n). 



s ,j, (n)- 



>G(z) 



■7> w (l| (n) 



> F(z) 



} 4-2 



“s (n) 



Figure 4.16— Block diagram of one-stage decomposition by a trous algoritlim. 

Now, the operation of decimation and then filtering s(n) resembles the portion of Mallat's 

algorithm, (4.53), where successively coarser approximations of s(n) were obtained by 
decimating and filtering each preceding approximation. For the signal detail is extracted merely 
by filtering the upper half of the spectrum below the Nyquist frequency. If the notation ^"(n) is 
adopted for the i'Mevel approximation of s(n) and w^n) for w(i, n), by induction, from (4.82) it 
follows that 

s (i+u (n) = f* s (i) (n) (4 g^ 

w (1) (n) = g * s (,) (n) 
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The block diagram for (4.83) is shown in Figure 4. 16. 



Earlier wavelet transform techniques resemble Gabor transforms in which sequences 
were decomposed by projection on representation vectors consisting of modulated Gaussian 
windows. Consequently, the wavelet function used by Shensa in his study of wavelets was a 
modulated Gaussian, also know-n as a Morlet window: 

(4.84a) 

Fourier transformation of (4.84a) produces 

\jt(co) = ^ JJk e -(co_v,: 2132 . (4.84b) 

In order to specify a wavelet function, it is, therefore, necessary to specify its center frequency v 

and its window rolloff factor p. 



Furthermore, in order to increase resolution in the Fourier domain it is possible to 

employ a wavelet function consisting of superimposed Gaussian envelopes modulated at 

different frequencies. In this case, the wavelet y(t) and its Fourier transform become 

M-l JSi(_!_y j v — ! — 

\|/(t) = I e ! ' - LM ) e - k,M 

k=0 

V|/(co) = | JTk I 2«e & ;l “ sl 

r U— 4 "i 



(4.85) 



Multiple frequency-domain translations of the Gaussian window are referred to as voices. 
Through the use of voices, the upper half of the frequency spectrum below the Nyquist 
frequency is partitioned into multiple spectral bins. Shensa listed a series of guidelines for 
designing a wavelet function with multiple voices: 
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The center frequency must lie in the upper half of the frequency spectrum below the 
Nyquist frequency. Consequently, 

7t/2 < v. (4.86) 

In order for \j/(t) to be an analytic function, 

(4.87) 

Obviously, (4.84b) has an infinite region of support. However, if, for (4.87), the equality 
is selected, 

y(0)= Jin e -2 ^ = 6.71 x 10 -9 . 

Consequently, the spectral content in regions of negative frequency will be negligible. 

To exclude aliasing, 

v<*-/2p. (4.88) 

The quantity Jl (3 is equal to one half of the Gaussian window where, at the edge of the 

passband, the frequency response of the filter is less than its maximum by a factor of e' 1 . 

Condition (4.88) ensures that the high-frequency edge of the passband is located below 

the Nyquist frequency. 

The fourth condition for designing wavelets employs the concept of relative bandwidth. 
The relative bandwidth B* of the window is defined to be the ratio of the window 
bandwidth to its center frequency, or 




Since, from (4.86) and (4.88), k/2 < v < iz, it follows that 

i^p = (3<B R <2P = ^. (4.89) 
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This condition is essentially a combination of previous conditions. 



5. Furthermore, if multiple voices are desired, it follows that the number of voices M 
required to completely cover the upper half of the spectrum is 

2 , 2(3 2)3 

As a result, an approximate choice for [3 is 

P = (4.90) 

6. To construct the separate voices, the mother wavelet \j/(t) is scaled by various factors 2 k/M 
where M indicates the number voices. Specifically, the k* voice v k (t) for k = 0, .... M-l, 

is given by 

*k«) = v(73r) (4.91a) 

from which it follows by the scaling property of the Fourier transform, 

v k (co) = 12™ JTn e-tf™®-") 2 '** 2 . (4.91b) 

This scaling produces the combined affects of increasing the Gaussian window decay 

factor by a multiple of 2 2 k/M and of translating the location of the window's maximum 

value to co = 2’ 2 k/M v. Selection of this scheme also maintains the affine nature of the 

transform. 

7. Then, to select the filter length, it is appropriate to consider the decay of the Gaussian 
window in the time domain. The filter will consist of a discretely sampled version of the 

continuous wavelet function. The DFT of the wavelet is, by definition, 

V k (e'“) = I v k (k)e-j“ m . 
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Due to the rapid decay of the window function V|/( t), the series can be truncated for values 



of the summation index k such that v M .,(k), the voice with the slowest time-domain rate 
of decay, will be negligible. If a negligibility threshold T u is specified, then a filter 
half-length K must be selected such that 

k„(K)| <T„. 



Inverting the envelope portion of (4.91a) produces for the filter half-length K 

k >± 2 < m -» m y-2 ln(T„) . 



(4,92) 



Consequently, the analyzing wavelet filter g(n) from (4.81) is defined to be 



g(n) = 



Vj/* (— n ) V n = -K, .... K 
0 otherwise 



(4.93a) 



and the filter g k (n) corresponding to the k' h voice v k (t) becomes 



gk(n) 



Vj/*(— n/2 k - M ) V n = -K, .... K 
0 otherwise 




(4.93b) 



Figure 4. 1 7-Spectral partitioning perfonned by three-stage decomposition by a trous method using four-voice 
Morlet window as analyzing wavelet. The interpolation filter, plotted as the lowpass filter represents a 
Lagrangian filter calculated from Daubechies' orthononnal scaling function on [0, 3]. 
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To illustrate the spectral behavior of the a trous method, Figure 4.17 plots the frequency 
decomposition performed by a three-stage structure employing a four-voice Morlet window. 

The spectral bin density increases in a logarithmic manner as digital frequency approaches zero 
and the width of each bin decreases similarly. The relatively large bandwidth of the lowest 
frequency bin occurs because of the relatively small number of stages. Additional stages would 
result in the continued trend of logarithmically increasing spectral bin density and bin width. 

The logarithmic trend must, however, be ultimately broken by the final spectral bin which will 
be defined by a half-band, lowpass filter. 

Figures 4.18a and 4.18b demonstrate the performance of the a trous algorithm for 
identification of the spectral components comprising a sequence. Figure 4.18a presents the 
decomposition of the lowpass test sequence (4.74a). Examining Figure 4.18a, the three 
harmonic components can be observed. At a scale of approximately 2.5 to 3, the band 
corresponding to the highest frequency component is evident. The middle frequency component 
is apparent in the range from three to four. The lowest frequency constituent is indicated by the 
sparce group of fins at scales of five to seven. The shape of the envelope is evident from the 
relative heights of the peaks of the two components possessing higher frequencies. However the 
lower frequency component appears relatively constant across the entire range of samples. This 
occurs because the integration time of the final scale extends the duration of the sequence. 

The a trous decomposition of the highpass test sequence (4.74b) is plotted in Figure 
4.18b. Unlike the Laplacian pyramid, the a trous method displays no difficulty in resolving the 
spectral components of (4.74b). Figure 4.18 provides an opportunity to understand the 
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relationship between the location of spectral components in the time-scale plane and their 
location in the frequency spectrum. The regions contained within each spectral bin are, again, 

indicated in Figure 4.17. The spectral component of the test sequence at 0.44-f s falls into the 
right-most spectral bin of Figure 4. 1 7. This component is indicated with in Figure 4. 1 8b at a 
scale of zero from approximately samples zero to 160. The location of the spectral component at 

0 . 2 1 - f s approximately coincides with the location of the boundary of the fifth and six spectral 
bins from the right-hand side of Figure 4.17. Since the analyzing wavelet contains four voices 



% 
K" ‘ 

% 

ro 




Figure 4.18a - Surface plot of a trous wavelet transform decomposition of 256-point, lowpass test sequence 
generated by (4.74a). Plot represents eight-stage decomposition using four-voice Morlet analyzing wavelet and 
Lagrangian interpolation filter calculated from filter for Daubechies’ orthonormal sealine function supported on 

[0, 3], 
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Figure 4.18b - Surface plot of a trous wavelet transform decomposition of 256-point, highpass test sequence 
generated by (4.74b). Plot represents eight-stage decomposition using four-voice Morlet analyzing wavelet and 
Laerangian interpolation filter calculated from filter for Daubechies' orthonormal scaling function supported in 

[0.3], 

(therefore there are four spectral bins per scale), the fifth and sixth spectral bins coincide with 
scales of 1.0 and 1.25, respectively. These scales coincide with an indication on 4.18b from 

samples zero to 96. Finally, the spectral component at 0.105-f. is located in the 10 th spectral bin 
of Figure 4. 1 7. This bin coincides with a scale of 2.25. Therefore, this final component, 
possessing the greatest power, provides an indication in Figure 4.18b throughout the entire range 
of samples. 

In working with the a trous algorithm, one significant disadvantage arises. Specifically, 
reconstruction of a sequence decomposed by the a trous method is extremely difficult. The 
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wavelet basis functions used in the a trous decomposition comprise frames in the sense described 
in Chapter II. Consequently, in order to reconstruct a sequence decomposed by the a trous 
method, it is necessary to construct biorthogonal set of bases. Daubechies outlined the details of 
construction the dual of a specified basis in [6], 

Apart from the difficult of inversion, the a trous algorithm, as implemented by Shensa, 
suffers from the disadvantage of long filter lengths. In the example under consideration, in order 
to satisfy a negligibility threshold (4.92) of 5 * 1 O' 3 , a filter length of 81 was necessary. Each of 
the four voices was implemented with a filter of equal length. Furthermore, the wavelet filter 
coefficients are complex. Complex filter coefficients increase the computational burden. 
Furthermore, from Figure 4.17 it is evident that the spectral decomposition resulting from the-a 
trous method is anything but power complementary. Partitioning in a manner similar to that of 
Figure 4.17 results in exaggeration of spectral content near the peaks of the voices and subdued 
indications between the peaks. In favor of the a trous algorithm, however, is the attainability of 
arbitrarily high resolution with smooth, regular filters. (It will be observed in the next section 
that filters sometimes suffer from lack of regularity.) 

E. MULTIRESOLUTION ALGORITHMS FROM CASCADED FILTER BANKS 

Section C of this chapter outlined the equivalence between Mallat's multiresolution 

analysis and a structure constructed from cascaded, two-channel, perfect reconstruction QMF 
banks. This section explores, in a more general manner, implementation of multiresolution 
analysis structures constructed from cascades of filter banks of arbitrary numbers of channels. 

For each method, the issues of division of the frequency spectrum and invertibility will be 
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considered. Spectral properties of the structures are of interest to obtain an understanding of the 
frequency resolution attainable with techniques under consideration. Invertibility is important 
because, in general, one good test of a decomposition technique is the accuracy of its 
reconstruction. 




a 3 (n) 



Figure 4, 79-Block diagram of three-stage, two-channel multiresolution structure constructed from cascaded QMF 
analysis banks. Algorithm represented is equivalent to Mallat's algorithm (4.53). 

A block diagram of the most basic multiresolution algorithm is depicted in Figure 4. 19. 
The process of Figure 4.19 is exactly equivalent to that of (4.53), Mallat's multiresolution 

decomposition. At each stage, the approximation sequence a k (n) is filtered with high-pass and 
low-pass half-band filters. The output of each filter is then dilated. The low-frequency channel 
decimator output— or the approximation channel, conforming to the terminology of 
multiresolution theory is then applied to another stage. The high-frequency, or detail, channel is 
employed for transmission, in the case of a communication system, or is applied to a detector in 
one possible case of a digital signal processing application. 



104 





If, hypothetically, the filter for each channel in the structure depicted in Figure 4.9 is an 
ideal filter, the DFT filter output sequence for the approximation sequence bandlimited to 

L 2 ([-7t/2,jc/ 2]), applying concepts from Chapter III. The DFT of the detail sequence is likewise 

bandlimited to L 2 ([7t/2, 3-7t/2]). Decimation of each sequence dilates its spectrum so that the 

approximation and detail sequence spectra are members of L 2 ([-7t, 7t]) and L 2 ([tc, 2-tc]), 
respectively. Consequently, the time-domain decimation and resultant frequency-domain 
dilation result, at each stage, in the examination of a smaller segment of the frequency spectrum. 




Figure 4-20--Partitioning of the frequency spectrum resulting from five-stage decomposition structure of Figure 
4. 19. Structure was implemented using, for the approximation and detail channels, the filters corresponding to 
Daubechies' orthonormal scaling function and wavelet, respectively, supported on [0, 13]. 

Figure 4.20 plots the spectral partitions resulting from a five-stage decomposition 
structured as in Figure 4.19. The spectral divisions indicated represent those obtained using the 
filters for Daubechies' orthonormal scaling function and wavelet on [0, 13]. The logarithmic 
contraction of succeeding spectral bins constitute the multiresolution analysis. The spectral 
partitioning of Figure 4.20 compares favorably to that obtained by the Laplacian pyramid in 
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Figure 4. 10. However, the localization provided by the partitioning of the frequency spectrum 
plotted in Figure 4.17 obtained by the a trous is obviously vastly superior to that of 4.20. For 
some applications, the resolution of Figure 4.20 may be inadequate while the computational 
burden necessary to obtain the results of Figure 4.17 is excessive. 

Brooks [33], in experimentation with signal processing techniques for ultra-wideband 
radar systems, observed that a dichotomous spectral decomposition, as is performed by the 
structure of Figure 4.19, provided poorer results than the a trous algorithm. One possible reason 
for the disparity between the performance of the two techniques considered was the inadequacy 
of the spectral resolution provided by Mallat's algorithm. Consequently, consideration of a way 
in which to improve the spectral resolution of multiresolution techniques without incurring the 
disadvantages of the a trous algorithm is justified. 








Figure 4.27~Multiresolution structure obtained by cascading detail channel of Figure 4.20 with a three-channel 

QMF-type filter bank. 

Figure 4.21 represents one possible approach to improving the results of Figure 4.20. 
Applying, at each stage, the detail channel to the input of the analysis bank of another filter bank 
produces a signal decomposition method capable of obtaining arbitrary precision. The study of 
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filter bank design is an advanced field and techniques for the design of perfect reconstruction 
filter banks of an arbitrary number of channels are available employing lattice structures [1 1] or 
through cosine modulation of a properly designed prototype filter [34], 

Figure 4.21 presents one example of a single stage of a structure using subchannel 
decomposition of the detail channel at each stage. In Figure 4.21, at each stage, the spectral bin 
corresponding to the detail sequence is further subdivided into three subchannels. Figure 4.22 
illustrates the spectral partitioning which results from a four-stage structure constructed from 
stages indicated by Figure 4.21. The improvement in spectral resolution is obvious. The density 
of groups of subchannels increases as digital frequency approaches zero. However, the overall 
structure does not provide the logarithmic spacing of spectral bins observed for the a trous 
decomposition as indicated in Figure 4.17. 




Figure 4. 2’-- Partitioning of frequency spectrum resulting from four-stage multiresolution structure constructed 
from stages appearing in Figure 4.22. Structure was implemented using, for the approximation and detail 
channels, the filters corresponding to Daubechies’ orthonormal scaling function and wavelet, respectively, 
supported on [0, 13]. Subchannel filter banks were implemented with three-channel filter bank coefficients are 

presented in [1 1]. 
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Subchannel decomposition can be performed using either perfect reconstruction QMF 
banks or pseudo-QMF banks— filter banks which do not strictly satisfy the perfect reconstruction 
criteria. Similarly, to the extent in which the pseudo-QMF bank deviates from the power 
complementary property, accurate representation in scale will be forfeited. 
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Figure 2. 2 j- Structure constructed from cascade of three-channel filter banks. 

An additional method to improve spectral resolution of multiresolution techniques was 

suggested by Gopinath and Burrus [35], Employing the term multiplicity M wavelet transform, 
Gopinath and Barrus suggested the concept illustrated in Figure 4.23 as a generalization of the 
dichotomous structure proposed by Mallat. In the frequency domain, spectral resolution 
increases as a base three logarithm instead of the base two logarithmic contraction observed in 
frequency spectra resulting from Mallat-type structures. Furthermore, structures similar to those 
of Figure 4.3 can be cascaded from filter banks of any number of channels. 
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Figure 4.24~Spectral partitioning resulting from four-stage decomposition constructed of structures as in Figure 
4.23. The filters indicated are 14-point pseudo-QMF filters designed in [l 1]. 

Figure 4.24 illustrates the partitioning of the frequency spectrum which results from a 
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Figure ‘/.25--Illustration of equivalence of two-stage decomposition-reconstruction structure with one-stage 

structure with delays in each channel. 

In reconstructing sequences decomposed by cascaded filter banks, the issue of delay must 
be addressed. As was discussed within the context of multirate system theory in Chapter II, a 



perfect reconstruction filter bank possesses an equivalent system transfer function which is equal 
to a delay. In decomposition-reconstruction structure "a" of Figure 4.25, the detail channel 



contains a delay L c representing, for example, a reconstruction filter bank used, as in Figure 
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4.21, to decompose the detail sequence into subchannels. The lower branch, the approximation 
channel, contains a two-channel bank corresponding to the next stage of decomposition. In the 
equivalent structure in structure "b" of Figure 4.25, the follow-on stage filter bank has been 
converted to the delay L c with which it is equivalent if it is a perfect reconstruction filter bank. 
Obviously, the delay in each channel must be equal. Otherwise, the approximation and detail 
sequences will be shifted with respect to each other and accurate reconstruction will not occur. 
Consequently, although a subchannel decomposition filter bank is not indicated in the detail 
channel of the final decomposition stage in Figure 4.25, if one were installed, an artificial delay 
in the approximation channel would be necessary. 
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Figure 4.2d--llustration of equivalent structures resulting from transmission of delay occurring in approximation 
and detail channels of a filter bank mulstiresolution structure. 

Figure 4.26 illustrates the total equivalent delay of a QMF decomposition-reconstructure 
containing a delay in each channel. In each channel a delay L s (equal to the order of the filters 
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used for decomposition into subchannels) is indicated which could include combinations of 
delays of subchannel decomposition-reconstruction structures, delays from subsequent 
decomposition-reconstruction stages and artificial delays introduced to resynchronize channel 

sequences. Since, by (3.9), [z ' L, ] T2 = z ' 2 , the total effective delay after the expander in each 
channel will be twice the delay of the channel prior to the expander. This equivalence is 
indicated by the transition from structure "a" of Figure 4.26 to structure "b." Furthermore, since 
the application of FIR filter and application of a delay represent linear operations, the two 
operations are also commutative. Consequently, as indicated in structure "c" of Figure 4.26b, 
the total, effective delay of the channel after the decimator can be moved beyond the summation 
point. Since a sequence reconstructed by a perfect reconstruction QMF bank is simply a delayed 
version of the original sequence, the total equivalent delay of structure "a" in Figure 4.26 is 

simply the sum of L c (the order of the filters used for decomposition into primary channels), the 

delay of the decomposition-reconstruction QMF bank structure itself, and 2 L S , the contribution 
from all of the delays introduced into the decomposed channels. 

To obtain a characterization of the total delay introduced into each channel of a 
multiresolution analysis constructed from cascaded filter banks, consider again structure "a" in 
Figure 4.25 as it would be implemented to reconstruct a sequence decomposed by M stages of 

the structure depicted in Figure 4.21. First, if a delay of L s , resulting from a subchannel 
decomposition-reconstruction structure, is introduced into the detail channel at stage M, of the 
final stage of decomposition, then an artificial delay of L s must be introduced into the 
approximation channel of the same stage to resynchronize the detail and approximation channel 
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sequences. Consequently, the total delay introduced into the detail and approximation stages 
consists of the sum of two contributions. The first source includes the delay L s resulting from 
decomposition-reconstruction structure corresponding to the decomposition of the detail 
sequence into subchannels and its corresponding equal delay in the approximation channel. 
Secondly, a delay L c corresponding to the primary channel decomposition-reconstructure at 
stage M must be considered. 

Continuing to the next stage, recombining the approximation and detail sequences from 
stage M-l to form the stage M-2 approximation sequence adds a delay of L c to the doubled 
equivalent delay present into the approximation sequence of stage M-l. Asa result, the total 
delay present in the approximation sequence at stage M - 2 is equal to 4-L s + 2-L c + L c . It is 

therefore necessary to introduce a delay of 4L S + L c + L c . to the detail sequence of stage M-2 
to resynchronize that stages approximation and detail sequences. 
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Figure 4.27 - Block diagram indicating structure necessary to obtain full expansion by channel from a sequence 
decomposed by M stages of the structure depicted in Figure 4.23. The values for L, and L, are defined in (4.94a) 

and (4.94b), respectively. 
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To obtain the structure required to obtain a full expansion by channel of the decomposed 
sequence by M stages of the structure presented in Figure 4.21, inductively the line of reasoning 
begun in Figures 4.25 and 4.26 must be pursued. The consequences of this induction appear in 

Figure 4.27. First, some arbitrary stage m, the subchannel sequence d m ,k) (n) must first be 
expanded and then filtered by the synthesis filter I k (z) corresponding to its subchannel. 

Artificial delays from two sources must then be introduced in order to re-synchronize the 
resultant stage m detail sequence component with its corresponding approximation sequence. 

The first delay Lj compensates for the artificial delay of L c introduced into the approximation 
sequence at stage M and is evaluated as 

L, = (2 M * m - l) L s . (4.94a) 

The second artificial delay contribution results from the primary approximation channel 
decomposition-reconstruction performed at all subsequent stages m+1 through M and is 
represented by a geometric series: 

M-m- 1 

L 2 =L c - I 2 P . (4.94b) 

p=0 

Next, the detail sequence component must be expanded and passed once through the 
primary detail channel synthesis filter. To complete the re-expansion, the sequence is subjected 
to m-1 iterations of expansion followed by filtering with the primary approximation channel 
synthesis filter. The final stage approximation sequence is fully expanded— as indicated by the 
second structure of Figure 4.27— by M iterations of expansion followed by filtering with the 
primary approximation channel synthesis filter. Application of the steps described by Figure 
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4.27 provides for a decomposed sequence a representation analogous to that depicted in Figure 
4. 1 1 for the Laplacian pyramid. The reconstructed sequence is simply the sum of all of the 
expanded channels or, if applicable, subchannels. 
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Figure 4.28a - Plots of reconstruction error versus number of decomposition stages for various decomposition 
stages applied to 256-point, lowpass test sequence produced by (4.74a). 

The accuracy of sequence reconstruction is dependent on the number of stages of 
decomposition performed and on the quality of the filter bank. For several stages of 



decomposition, small errors resulting from filter imperfections, such as roundoff error, 



accumulate and corrupt the reconstruction process. To obtain optimum results, filter banks 



which strictly satisfy the perfect reconstruction criteria should be used. In using pseudo-QMF 



filter banks which do not satisfy the perfect reconstruction, poorer results will be obtained. 



These concepts are illustrated by Figures 4.28a and 4.28b. Figure 4.28a indicates reconstruction 



eiTors for multiresolution structures applied to the test sequence generated using (4.74a). 



The lowpass test sequence from (4.74a) was constructed from low-frequency harmonics. 



From the spectral partition diagrams presented for the various multiresolution structures 
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considered, the spectral resolution of the multiresolution structures was greatest in the lower 
regions of the frequency spectrum. In Figure 4.28b, the reconstruction error for the highpass test 
sequence (4.74b) is presented. Illustrating results for the test sequence generated by (4.74b), 
Figure 4.28b illustrates multiresolution technique performance for signals constructed from 
frequency components located in the upper half of the frequency spectrum below the Nyquist 
frequency. 
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Figure 4.28b - Plots of reconstruction error versus number of decomposition stages for various decomposition 
stages applied to 256-point test sequence produced by (4.74b). 

In the case of both test sequences, the lowest reconstruction errors resulted from 



multiresolution structures comprised strictly of perfect reconstruction filter banks. In 



two-channel, zero-subchannel and two-channel, two-subchannel structures, the filters used were 



all obtained from Daubechies' orthonormal wavelet and scaling function supported on [0, 13], 



Consequently, the multiresolution analyses obtained were equivalent to true, orthogonal 



decompositions. 
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For the two-channel, three-subchannel, the two-channel, four-subchannel and the 
three-channel, zero-subchannel cases, the structures all entailed the use of pseudo-QMF banks 
which did not strictly satisfy the perfect reconstruction criteria. The structures constructed using 
the three-channel filter banks employed the set of filters designed by Vaidyanathan in [11]. The 
four-channel filter bank used a set of 30-point, fourth-band filters designed by spectral 
factorization of filters obtained by the McClellan Parks algorithm. The design process for the 
four-channel filter bank will be outlined in Chapter V. Because of failure to satisfy the perfect 
reconstruction property, structures using these methods consistently produced greater 
reconstruction error. 

TABLE 4. 7-PARTITION OF SPECTRUM RESULTING FROM EIGHT- ST AGE, TWO-CHANNEL, • 
ZERO-SUBCHANNEL MULTIRESOLUTION DECOMPOSITION. 



Decomposition stage 


Scale 


Spectral bin 
lower bound 


Spectral bin 
Upper bound 


Stage 1 Detail 


0 


2 H 


2‘H 


Stage 2 Detail 


1 


2‘H 


2‘H 


Stage 3 Detail 


2 


2 'H 


2’H 


Stage 4 Detail 


3 


2‘H 


2 H 


Stage 5 Detail 


4 


2' 6 f s 


2 H 


Stage 6 Detail 


5 


2'H 


2‘ 6 -f s 


Stage 7 Detail 


6 


2‘H 


2 H 


Stage 8 Detail 


7 


2'H 


2' 8 -f s 


Stage 8 Approximation 


8 


0 


2 H 



Inductively extending the results of the spectral partition diagram (4.20) corresponding to 
this structure provides insight into the relationship between scale and frequency. The first 
spectral bin, which coincides with scale zero in Figure 4.29a, constrains the region of the 
frequency spectrum from 0 . 5 • f s to 0.25-f $ . A scale of unity coincides with the region of the 
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frequency spectrum from 0.25 f s to 0.125-f s , the region corresponding to the signal detail 
extracted during the first stage of decomposition.. Inductively, therefore, a scale of m, 
corresponding to the signal detail extracted at the m th stage of decomposition, coincides with the 
spectral bin partitioning the region from 2' <m ~ l) -f s to 2' (m * :) -f s . This subdivision is tabulated in 
Table 4.1. 
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Figure 4.29a — Surface plot of decomposition of lowpass test sequence (4.74a) using eight-stage, two-channel, 
zero-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthonormal wavelet and scaling function supported on [0, 13]. 

Figures 4.29a and 4.29b demonstrate the full expansion by channel of decomposition of 
sequences (4.74a) and (4.74b), respectively. In Figure 4.29a, the depiction of the lowpass 
sequence (4.74a), the resolution of scale of the decomposition at low-frequencies is apparent. 



The component of sequence (4.74a) located at digital frequency 2 -jt 9/512 (= 0.0 1 8-f s ) falls in 
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the range of 2°f s to 2' 6 -f s which corresponds to a scale of four. In Figure 4.29a, a series of 
circular contours centered on the scale axis for m = 4 can be observed from approximately 
samples 32 to 216. Similarly, the spectral component located at digital frequency 2-71-31/512 

(= 0.06-f s ) lies in the region 2‘ 5 -f s to 2' 4 f s corresponding to a scale of m = 3. Asa non-integer 

power of two, the location 0.06-f s = 2‘ 4 05 which is very near the boundary between adjacent 
scales two and three. Consequently, much of the component’s energy appears spread between 
scales two and three. The component at digital frequency tc/ 256 (=0.002-^ exists at the 
boundary between scales six and seven. The indication of this component, which completes 
only a single oscillation during the duration of the test sequence, appears spread between the two 
scales. 

Figure 4.29b presents the two-channel, zero-subchannel decomposition of highpass test 
sequence (4.74b). The relative coarseness of the multiresolution decomposition in the higher 
regions of the frequency spectrum becomes apparent in Figure 4.29b. Because of the rapid 
variations with respect to time, the contour plot in Figure 4.29b is difficult to read. However, 
information can be gained from examining the surface plot. Again, comparing the plot with the 
spectral partition diagram, the spectral component of (4.74b) which lies at a digital frequency of 

2-71-27/512 (= 0.105-fJ falls within the region of 2' 3 to 2' 4 corresponding to a scale of m = 2 as 
plotted in Figure 4.29b. The scale of m = 2 corresponds to the range of peaks with highest 
amplitude. The spectral component located at 2 -tt 55/256 (= 0.215* f s ) is contained within the 
region 2 * 2 to 2' 3 which corresponds to a scale of m = 1 as plotted in Figure 4.29b. The 
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indication of this component is recognizable in the second range of peaks between the highest at 
m = 2 and the range that runs along the edge of the surface plot. Finally, the spectral component 

at digital frequency 2 - k - 1 1 3/256 (~ 0.44-f s ) lies within the region of 2' 1 to 2 ' 2 and is indicated in 

the range of peaks at scale m = 0 which runs along the edge of the surface plot. For the highpass 
test sequence, each of the harmonic components can be resolved. However, in the case of the 

highest frequency components, the spectral separation is approximately 0.2-f s . Furthermore, 
each of the three spectral components fell within adjacent spectral bins. Consequently, the 
example demonstrates the maximum resolution for a two-channel, zero-subchannel 
multiresolution structure. 
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Figure 4.29b — Surface plot of decomposition of highpass test sequence (4.74b) usmg eight-stage, two-channel, 
zero-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthonormal wavelet and scaling function supported on [0, 13]. 
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Figure 4. 30u--Surface plot of decomposition of lowpass test sequence (4.74a) using eight-stage, two-channel, 
three-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubeckies' orthonormal wavelet and scaling function supported on [0, 13]. The subchannel decomposition was 

performed with pseudo-QMF filters designed in [1 1], 

Figures 4.30a and 4.30b provide two-channel, three-channel decompositions of the 
lowpass test sequence (4.74a) and the highpass test sequence (4.74b), respectively. From the 
spectral partition diagram for this structure. Figure 4.22, each spectral region 2' (nW) to 
is further subdivided into three subregions. For the general case, a multiresolution structure 
constructed to divide, at each stage, the detail sequence into M subchannels, at the m* stage of 
decomposition, the k* subchannel corresponds to the region of spectrum contained in the interval 
^2 -(m+2) + (2~( m+, )-2- <m+2) ) j • f s , ^2 -(m+2) + -jj-j- ■ (2 -(m+1) -2 -(m+2) ) j • f s . 

Because of the increase in the density of spectral subdivision afforded by decomposition 
of the lowpass test sequence into subchannels, the time-scale localization presented in Figure 
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4.30a provides time-scale localization which is improved over results presented in Figure 4.29a. 
In particular, the component at digital frequency 2-^-3 1/5 12 clearly localized at a scale of three 
instead of being partially spread over adjacent scales as occurred in Figure 4.29a. Similarly, the 
component at digital frequency 2-71-9/5 12, which roughly corresponds to a digital frequency of, 

c o 

2' -f s , can be recognized to exist primarily between scales of four and Five. The component 
characterized by the lowest frequency is still recognizable in the vicinity of scales six and seven. 
However, an indication appears between scales of five and six which cannot be accounted for 
based on the components known to exist in the test sequence. This unaccounted for component 
may have appeared as a result of spreading from the component at scale six. 
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Figure 4. 306- Surface plot of decomposition of highpass test sequence (4.74b) using eight-stage, two-channel, 
three-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies’ orthonormal wavelet and scaling function supported on [0, 13]. The subchannel decomposition was 

performed with pseudo-QMF filters designed in [1 1]. 
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In the case of two-channel, three-subchannel decomposition, the spectral components of 
the highpass test sequence (4.74b) have also become more distinct. In the plot of Figure 4.30b, 
because the spectral content of the test sequence has been localized in specific subchannels and 
excluded from others, resulting in an indication of three separate components. Consequently, the 
indication provided by the contour plot constitutes a more useful display than the contour plot 
for Figure 4.29a. 
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Figure 4. 3 1 a--Sur face plot of decomposition of lowpass test sequence (4.74a) using eight-stage, two-channel, 
four-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthonormal wavelet and scaling function supported on [0, 13]. The subchannel decomposition was 
perfonned with pseudo-QMF filters obtained from spectral factorization of fourth-band filters designed by the 

McClellan-Parks method. 

Examining Figure 4.31a, a plot of a two-channel, four-subchannel decomposition of 
lowpass test sequence (4.74a), increasing structural complexity from three to four subchannels 
provides some improvement in localization for the test sequence presented. The component of 
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digital frequency 2 rc-9/5 12 appearing between scales four and five have become distinct from 
the indication unaccounted for between scales five and six. Consequently, it is likely that the 
false indication resulted from spillover of energy from the component at 2 ^-2/512. Between 
scales five and six, the resolution of the decomposition using four subchannels becomes 



0.25 -2' 6 f s (=f/5 12). Beneath those scales the resolution continues to increase logarithmically in 
base two. It appears, therefore, as though the four-subchannel decomposition may provide 
excessive resolution for the test sequence demonstrated. In general, the principal benefits of 
decomposing sequence into subchannels occurs in the higher regions of the frequency spectrum. 
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Figure 4 J/ 6-Surface plot of decomposition of highpass test sequence (4.74b) using eight-stage, two-channel, 
four-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies’ orthonormal wavelet and scaling function supported on [0, 13]. The subchannel decomposition was 
performed with pseudo-QMF filters obtained from spectral factorization of fourth-band filters designed by the 

McClellan-Parks method. 
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Figure 4.3 lb further illustrates the principal afforded by subchannel decomposition. 
Specifically, in Figure 4.31b, an improvement in the resolution of the upper regions of the 
frequency spectrum is evident. In Figure 4.31b, a plot of the two-channel, four-subchannel 
decomposition of highpass test sequence (4.74b), the sequence component at digital frequency 
2-tt 1 13/256 has been localized as distinct from scale zero. Furthermore, as compared with the 
two-channel, three-subchannel plot of figure 4.30b, the test sequence component at digital 

frequency 2 -jt 27/256 exhibits improved localization. 

TABLE 4. 2-PARTITION OF SPECTRUM RESULTING FROM FIVE-STAGE, THREE-CHANNEL, 
ZERO-SUBCHANNEL MULTIRESOLUTION DECOMPOSITION. 



Decomposition 

stage 


Multiresolution 

scale 


Spectral bin 
lower bound 


Spectral bin 
upper bound 


Stage 1 


Scale 0.0 


f s /3 






Scale 0.5 


f s /6 


f s /3 


Stage 2 


Scale 1.0 


f s /9 


f s /6 




Scale 1.5 


f s /18 


f,/9 


Stage 3 


Scale 2.0 


f s /27 


f/18 




Scale 2.5 


f s /54 


f s /27 


Stage 4 


Scale 3.0 


f s /81 


f s /54 




Scale 3.5 


f/162 


f s /81 


Stage 5 


Scale 4.0 


f/243 


f/162 




Scale 4.5 


f/486 


f s /243 




Scale 5.0 


0 


f s /486 



As indicated previously, decomposition into three primary channels provides base three 
logarithmic partition of the frequency spectrum. Specifically, the first stage of decomposition 
divides the frequency spectrum into spectral bins covering the regions f $ /6 to 2 f s /6 and 2-f./6 to 
3-f,/6. During the second stage of decomposition, the lower third of the frequency spectrum is 
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partitioned along intervals f/18 to 2 • f s / 1 8 and from 2 • f s / 1 8 to 3 f s /18. At stage m, the signal 
detail extracted will occupy the spectral regions f s /(2-3 m ) to 2-f s /(2-3 m ) and from 2-f s /(2-3 m ) to 

3-f s /(2-3 m ). Table 4.2 presents tabulated information regarding the spectral subdivision 
occurring from a three-channel, zero-subchannel multiresolution decomposition. 
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Figure 4J2a--Surface plot of decomposition of sequence (4.74a) using five-stage, three-channel, zero-subchannel 
multiresolution structure constructed from cascaded QMF banks. The decomposition was perfonned with 

pseudo-QMF filters designed in [1 1], 

Finally, Figures 4.32a and 4.32b present the results of three-channel, zero-subchannel 
decomposition of the lowpass test sequence (4.74a) and highpass test sequence (4.74b), 
respectively. Because of the base three logarithmic partitioning of the frequency spectrum, 
localization of sequence spectral components in the lower portion of the frequency spectrum will 
be improved over the localization resulting from a two-channel decomposition. Specifically, for 
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sequence (4.74a), component at digital frequency 2-71-2/512 is located, according to Table 4.2, in 
the vicinity of the boundary between scales 4.0 and 4.5. In fact, most of the energy of this 
component appears in the channel corresponding to scale 4.0. Similarly, the component at 
2 • 7T - 9/5 12 is located in the vicinity of the boundary between scales 2.5 and 3. The energy of the 
appears mostly at scale 2.5. Finally, the energy of the component at digital frequency 
2-7T-3 1/512 is indicated at scale 1 .5. An unexplainable indication appears at a scale of 0.5 
between samples zero and 32 and also at approximately sample 96. From Figure 4.24, the filters 
from which the three-channel filter bank was constructed produced greater stopband error than 
the two-channel filter bank. Consequently, energy from the sequence appears to leak, producing 
less precise localization of components of the test sequence. 

With regard to Figure 4.32b which presents the three-channel, zero-subchannel of the 
highpass test sequence (4.74b), the structure appears to localize sequence components 
approximately as well as the two-channel, zero subchannel. Comparing the two-channel, 
zero-subchannel frequency partition diagram Figure 4.20 with the corresponding diagram Figure 
4.24 for the three-channel, zero-subchannel structure, the two methods provide very similar 
resolution in the upper regions of the frequency spectrum. Specifically, the three-channel 

structure divides the spectral region between f s /2 and £76 into two channels while the 
two-channel structure divides the region between f s /2 and fj/4 into two channels. This similarity 
is reflected in the Figure 4.32b. The component at digital frequency 2-7t-27/256 appears as the 
largest range of peaks at a scale of two. The component at digital frequency 2-tt-I 13/256 is 
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indicated at a scale of zero along the edge of the surface. The indication of the component at 
2-3X-55/256, however, is less well localized. Indications of this final component appear to be 
present at scales of 0.5 and 1.0. 
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Figure 4.32b--Suiface plot of decomposition of highpass test sequence (4.74b) using five-stage, three-channel, 
zero-subchannel multiresolution structure constructed from cascaded QMF banks. The decomposition was 

performed with pseudo-QMF filters designed in [1 1], 

Finally, the three-channel structure demonstrated in Figure 4.32a presents a noteworthy 
phenomenon. The surface in the higher scale (lower frequency) regions possesses a jagged 
texture. When the surface resulting from expansion by a filter is characterized by rapid 
variations such as are shown in Figure 4.32a, the filter is said to lack regularity. The appearance 
of this jagged quality is noteworthy because it occurs in the high-scale channels containing the 
slowly varying components. For the the two-channel structures, the surfaces in the regions of 
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higher scales presented relatively smooth appearances. The non-regular surface observed in 
Figure 4.32a increases the difficult in reading of the surface plot. The issue of regularity will be 
addressed again in Chapter V. 

To provide an overall assessment of the three-channel structure, it affords greater 
resolution at higher scales (lower frequencies) without incurring the increase in structure 
complexity involved in performing decomposition into subchannels. Compared to the 
two-channel, zero-subchannel structure, the structure provided somewhat improved resolution. 
However, the performance in the high-frequency regions is not as great. Consequently, included 
among the issues in deciding whether to employ three channels or two-channels with 
subchannels is whether resolution at higher frequencies is critical. 
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V. BASIS FUNCTIONS AND FILTERS FOR .MULTIRESOLUTION 

DECOMPOSITION 

A. INTRODUCTION 

Chapter IV developed the theory of multiresolution decomposition and demonstrated its 
equivalence to the cascading of QMF bank structures described in Chapter III. In Chapter IV, 
the wavelets and scaling functions used for multiresolution signal decomposition were addressed 
only the most abstract sense. In the present chapter the nature of these functions will be 
addressed in greater detail. Section B outlines some of the basic concepts of the theory of 
two-scale difference equations, the class of equations from which wavelets and scaling functions 
were, in Chapter IV, shown to be derived. In Section C, some of the properties of Daubechies' 
orthonormal wavelets and scaling functions will be considered. Finally, Section D will return to 
a filter bank perspective for the study of multiresolution decompositions. Some basics of the 
theory of perfect reconstruction and pseudo-QMF bank filters will be addressed. 

It must be noted that a number of approaches exist for the construction of wavelets and 
scaling functions which will not be addressed in the present study. One of the earliest methods 
w hich has subsequently become relatively well-developed is to construct wavelet and scaling 
function filters from polynomial splines [36, 37], In the construction of w'avelets from splines, a 
polynomial spline is used for the scaling function. Since polynomial splines of non-zero order 
are not orthogonal with respect to integer translation, the scaling function family constitutes a 
frame in the sense described in Chapter 2. The biorthogonal basis functions which complement 
the polynomial splines to complete the frame can be constructed from, among other methods, 
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using recursive (HR) filtering techniques. Finally, the design of wavelets has been approached 
from an optimization theory perspective, as well [38]. 



B. TWO - SCALE DIFFERENCE EQUATIONS 

In equation (4.23) (repeated here for convenience), it was shown that scaling functions 

represent the solutions to two-scale difference equations: 

■f • <Ht) = ^ h(k) • 0(u - k) . (5.1) 

k 

Before attempting to solve (5.1), it is first useful to consider some of its properties. The 
development which follows makes the assumption that (f> is a real function. However, the results 
can easily be generalized to complex functions as well. First, integrating each side of (5.1) over 
its region of support produces [27] 

■r • J<Mt) du = 1 h(k) • j 0(u — k) du. (5.2) 

k 

Applying the changes of variables v = u/2 and w = u -k to the appropriate sides of (5.2) results in 

J 4>( v) dv = Z h(k) • } (J)(w) dw. (5.3) 

k 

Now, each side of (5.3) contains a common factor of J <f>(u) du, which, if it is non - zero , can be 

cancelled on both sides of the equation resulting in 

I h(k) = 1 . (5.4) 

k 

Two properties have now been specified for the families of scaling functions which 
represent solutions to two-scale equations. The first is (5.4) which indicates that the filter whose 
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impulse response is h( k) possess a DC gain of unity. The second, which was necessary to obtain 
(5.4) requires that 

\ <j>(u) du * 0. (5.5) 

In fact, if the objective is to construct orthonormal basis functions, it is desirable to specify that 

1 0(u) du = 1. 

If the condition is imposed that the family {0(u - k)} kg 2 constitutes an orthonormal 

basis, then another interesting property for h(k) can be obtained. First, the change of variables t 
= 2-u is made for (5.1) resulting in 

0(t) = 2 • Xh(k) • 0(2 • t-k). (5.6) 

k 

Now, if 0 is orthonormal with respect to integer translation, then j 0(t -k) 0(t - j) dt = 5j k where 

5 is Kronecker's delta function. Furthermore, j q>( 2-t -k) <J)(2 t - j) dt = k 2 also holds. 

Consequently, evaluation of the projection of integer translates of 0 onto each other yields 

j 0(t-m) • 0(t-n) dt = 4 ■ XI h(k) • h(j) • J 0(2 • (t - m) -j) • 0(2 • (t- n) -k) dt. (5.7) 

k j 

Applying the change of variables u = 2-(t - m) to (5.7) produces 

J 0(t — m) • 0(t - n) dt = 2 • LX h(k) • h(j) • J 0(2 • u -j) • 0(u + 2 • (m- n) -k) dt. (5.8) 

k j 

Next, applying the consequences of the assumption of orthonormality to both sides of (5.8), 
results in 

6 m . n =2-ZZh(k)-h(j)-Sj.k-2<m-n,. (5-9) 

k j 

Finally, by the summation sifting property of the Kronecker delta function, 
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(5.10) 



5 m .n = 2 • X h(k) • h( k - 2 • (m- n)). 

J 

Equation (5.10) provides the basis for at least two conclusions. First, for the case when 

m = n, 

X h 2 (k) = y . (5.11) 

k 

In other words, if the two-scale difference equation solution § is orthonormal with respect to 
integer translation, then the sum of the squares of the coefficients of the corresponding 
generating filter is one half. 

A second interpretation, made by Vetterli [25], state that the impulse response h(k) of the 
generating filter is orthogonal respect to even translates of itself. This is equivalent to another 
interpretation. Define the sequence g(n) from the convolution of h(n) with its own time 
reversal h(n): 

g(n) = 2 • X h(k) • h(n-k). (5.12) 

k 

If g(n) is defined in this manner, then (5.10) is equivalent to the expression 

g i2 (n) = g(2n) = 5 0 , n . (5.13) 

Now, (5.13) is equivalent to (4.78). Therefore, the filter g(n) is an a trous filter in the sense 
described by Shensa [31] or a half-band filter in the sense described by Mintzer [32], 
Furthermore, (5.13) requires that the length of h(k) be an even integer [25]. If the length of h 
were an odd integer L, then g(L - 1) = 2 h(0)-h(L-l) *0, and (5.13) would not hold. Finally, the 
evaluating the Z-transform of (5.12) results in 

G(z) = H(z) • H(z) = H(z) ■ H(l/z). (5.14) 
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By (5. 14), H(z) is simply a spectral factor of a half-band filter. This concept will be employed 
later in Section D of this chapter within the context of designing filters for QMF banks. 

Solving for scaling equations given an appropriate filter can be performed in either the 
time domain or in the Fourier domain. The Fourier-domain solution begins with Fourier 
transformation of (5.1), which was shown in (4.36b) to produce an expression equivalent to 

p((o) = H(e J f).0(^). (5.15) 

By induction, (5.15) produces [27] 

p(co) = p(0)- ft Ffiej®' 2 *). (5.16a) 

k= 1 

Now, if the family of scaling functions p is orthonormal, then j p(u) du = 1 which is equivalent 

A • 

to p(0) = 1 . Therefore, (5.16a) becomes 

p(d))= flH(eJ 01/211 ). (5.16b) 

k= 1 

The Fourier transform of the scaling function p is simply the infinite product of progressively 
more dilated versions of the Discrete Fourier Transform (DFT) of the generating filter h(k). 
Moreover, as the product index k tends towards infinity, the DFT term H(e j ^ ) on the 
right-hand side of (5.16b) becomes infinitely spread out in the Fourier domain. The Fourier 
transform of p, therefore, has an infinite region of support which indicates the possibility of 
finite support in the time domain. Compact support of the scaling function p is highly desirable 
since, because of how h(k) was defined in equation (4.23), only a compactly supported scaling 
function will produce a finite impulse response filter. 
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To solve the scaling function in the time domain, it is first convenient to absorb the 
factor of two on the right hand side of (5.6) into the filter coefficients by defining h {2) (k) = 
2h(k). Consequently, an equivalent form for (5.6) becomes 



9(0 = Ih, 2 )(k) • 9(2- t-k). (5.17a) 

k 

The time-domain solution of (5. 1 7) begins by solving for values of the function <J> at integer 
points in its domain [39]. Consequently, the equation to be solved becomes 



9(n) = Zh ( 2)(k) • <5(2 • n - k). 
k 



(5.17b) 



Now, if the scaling function <J> is strictly supported on the region [0, 2 L - 1] and its generating 
filter h(k) is of length 2-L - 1 for some inter L, then (5.17b) can be expressed in matrix form as 






<t>(0) 

9(D 

<5(2- L-l) 



\ 



h ( 2)(0) 0 0 

h<2)(2) h(2) ( 1 ) h(2)(0) 



0 

0 



0 

0 



h<2)(2 • L - 1) h <2) (2 • L - 2) h(2)(2 • L - 3) 

0 0 h (2) (2 L-l) ) 



m 

<5(D 

9(2- L-D 

(5.18) 



It is obvious that (5.18) is simply an eigenvalue problem. The vector containing the values of 
the scaling function 9 for integer arguments is simply the eigenvector of a unit eigenvalue for the 
double-right-shift matrix of filter coefficients in (5.18). 



Now that the value of the scaling function 9 has been solved for integer arguments, the 
next step is to rewrite (5.17b) in a form that will yield half-integer values. 

9(|) = Ih, 2) (k)-9(n-k). (5.18) 

z k 
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Equation (5. 18) is obviously a simple convolution, nevertheless, insight is gained from 
expressing it in matrix form: 



o \ 

0(t) 



J 



' <t>(0) 


0 


o N 


0(1) 


0 


0 


0(2) 


0(0) 




0(2 • L - 1) 


0(2 • L - 3) 


0(0) 


0 


0(2 -L -2) 


0(1) 


, 0 


0 


0(2 - L- 1) J 



h«2)(0) 

h (2 )(0) 

h { 2)(2 * L- 1) 



(5.19) 



From (5.19), a pattern begins to emerge. Extended (5.1 8) to the general dyadic case produces 

■K^) = Ih 0 )(k)-<i>(£-iO. (5.20) 

Finally, adopting the vector notation 

A(ra)_f JL _!_ ... 

t 1 2 m 2 m 2 m 



h= h (2) (0) h, 2) (l) ••• h ( 2)(2 • L - 1) 



and 






0(£) 


0 


0 > 


0(£) 


0 


0 


Wr$r) 


<K^r) 


o 


0( 1 S tI ) 




o 


A/ (2L-l)-2 m x 
0( 2 m ) 


j./ (2 L-2)-2 m , 
0( 2 ™ ) 


0(£) 


0 


x, (2 L-2)-2 m +l , 
0( 2 ra ) 


0(£) 


, 0 


0 


••• 0( 2 m ) J 
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equation (5.20) can be expressed in a more compact form: 



^ m *i) = <D<nu h. (5.21a) 

Similarly, if the vector \j/ m) for the wavelet function y(t) is defined in a manner analogous to 
<J> (inl for the scaling function 4>(t) and the wavelet generating filter g in in a manner analogous to 
the scaling function generating filter h, then, as a consequence of (4.51c), the values of \|/(t) can 
be obtained at dyadic points from 

\|/ m+1) = <D tm) g. (5.21b) 




Figure 5./a--Function generated by four-stage dyadic expansion of a 20-point lowpass filter with cutoff frequency 
of 0.375 f s . The function represents an example of a function which lacks regularity. 

The recursive operations described by (5.21a) and (5.21b) represent dyadic expansions of 
the functions <f)(t) and \y(t). At each iteration the discrete sequences characterized by the vectors 

and \|/ m> are expanded by a factor of two and then convolved with filters h(n) and g(n), 
respectively, whose coefficients are contained in h and g. After m iterations, the values of the 
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of the scaling function $(t) and wavelet \|/(t) will be known at all points t n in region of support 
such that t n = n/2 m . 

A final issue which must be addressed within the scope of two-scale difference equations 
is the topic of regularity. Regularity is defined to be the degree to which a recursively expanded 
function, such as a two-scale difference equation, converges to a continuous smooth function 
[26]. The conditions under which the regularity of a function can be guaranteed of themselves 
constitute grounds for an extensive study. However, one apparent example of a requirement for 
regularity was observed during experimentation performed during the course of the present 
study. For regularity, it is necessary that, for the filter to be expanded, the DFT exist primarily 
within the lower half of the frequency spectrum below the Nyquist frequency. Figures 5.1a and 
5.1b show two examples of the dyadic expansion of two candidate filters for the generation of a 
scaling function. Both filters were designed by performing spectral factorization of lowpass 
filters designed using the window method [12]. For the first filter the cutoff frequency was 

selected as 0.375-f s and the second filter was assigned a cutoff frequency of 0.25-f s . As can be 
seen in Figure 5.1a, the function possesses an erratic appearance. As the bandwidth of the 
candidate function generating filter increases, the generated function becomes less smooth and 
more and more erratic. 

On the other hand, Figure 5.1b shows a function generated through dyadic expansion of a 
filter designed with exactly the same specifications except that the cutoff frequency was 
specified as 0 . 2 5 * f s . The plot of the function generated by the second filter displays significantly 
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greater smoothness. Figure 5.1b represents an example of a function which, as far as can be 
determined from the plot, possesses the quality of regularity. 




Figure 5 J 6— Function generated by four-stage dyadic expansion of a 20-point lowpass filter with cutoff frequency 
of 0.5 f s . The function represents an example of a function which possesses regularity. 

C. DAUBECHIES * FAMILY OF ORTHONORMAL WAVELETS 

Ingrid Daubechies of Bell Laboratories made a landmark contribution to the field of 
multiresolution analysis when, in 1988, she published a mathematical paper on orthonormal 
wavelets [21]. Aside from providing a general overview of the progress which during recent 
years had occurred within the field of wavelets, Daubechies' paper described the construction of 
a class of functions which possessed many of the most desirable aspects of wavelet functions. 

The functions were specified in terms of their generating filters. Even-length filters from length 
four to 20 were tabulated in the paper. Each of the filters satisfied (5.4) and (5. 1 1). 

Furthermore, they each maximized the number of zeros on the unit circle at z = -1, a condition 
which ensures maximum regularity [26], In this section plots of two wavelet-scaling function 
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(pit), Ip < t ) 



pairs will be presented, and one will be examined more thoroughly from the perspective of 
implementation of the multiresolution structures from Chapter IV. 




Figure 5.2a--Superimposed plots of Daubechies' orthonormal wavelet and scaling function on [0, 3]. Plots were 
generated by four-stage dyadic expansion of corresponding filter listed in [21]. 
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Figure 5. 2b— Polar plot of zero locations for generating filter for Daubechies* orthonormal wavelet and scaling 

function supported on [0, 3]. 
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Figure 5.3— Two-channel, two-subchannel multiresolution decomposition of 256-point lowpass test sequence 
(4.74a). Decomposition was performed using generating filter for Daubechies' orthonormal basis function 

supported on [0,3]. 

Perhaps the most famous of Daubechies' orthonormal scaling function-wavelet pair is 
plotted in Figure 5.2a. The functions are generated from a four-point filter and are supported on 
[0, 3], As can be inferred from viewing the plots, the smoothness of these functions is marginal. 
Figure 5.2b provides a polar plot of the zeros of the generating filter. Of its three poles, two fall 
on the unit circle. Finally, as indicated in Chapter IV, the wavelet transform, which essentially 
consists of the detail function (4.33), is comprised of a plot of weighted, shifted versions of the 
wavelet. Consequently, it is expected that the texture of a surface plot of a wavelet 
decomposition should reflect the shape of the wavelet from which it is constructed. Figure 5.3 
supports this expectation. In particular, the peaks in the which appear at higher scales reflect the 
shape of the wavelet function plotted in Figure 5.2. 
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Figure 5.4 — Superimposed plots of Daubechies' orthonormal wavelet and scaling function on [0, 13]. Plots were 
generated by four- stage dyadic expansion of corresponding filter listed in [21]. 




Figure 5. 5— Polar plot of zero locations for generating filter for Daubechies' orthonormal wavelet and scaling 

function supported on [0, 13]. 

Figure 5.4 presents superimposed plots of Daubechies' orthonormal basis function and 
scaling function supported on [0, 13]. This pair was used extensively throughout Section IV. E 
to illustrate the performance of multiresolution decomposition techniques. The zeros of the 
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generating filter are plotted in Figure 5.5. Of the 13 zeros of the generating filter’s 
characteristic polynomial, seven occur at the location z = -1. 




L4E-012 

1.2E-012 

IE-012 

8E-013 

EE-013 

4E-013 

2E-013 

0 

-2E-013 

-4E-013 

-EE-013 



i 

cr 



I 



\ 



Figure 5. 6— Time-domain demonstration of performance of QMF bank constructed from generating filter for 
Daubechies' orthonormal wavelet and scaling function on [0, 13]. Left-hand axis indicates equivalent impulse 
response while right-hand axis indicates deviation t(n) - 8(n - N slax ) from ideal perfect reconstruction 

performance. 

In Section IV.E, multiresolution decomposition techniques were addressed from the filter 
bank perspective. In Chapter III, the degree to which specific examples of filter banks satisfied 
the perfect reconstruction criterion and the power complementary property was demonstrated. 
The question, therefore, arises as to the degree to which the generating filters for Daubechies' 
orthonormal wavelets and scaling functions satisfy the perfect reconstruction and power 
complementary properties. Figures 5.6 and 5.7 illustrate the performance of these functions 
employing the same plots used in Chapter III. Figure 5.6 presents a plot analogous to Figure 
3.14. Upon comparing Figure 5.7 with Figure 3.14, it is apparent that Daubechies' generating 
function provides vastly superior performance for filter bank applications. The filter bank 

demonstrated in Figure 3.14 deviated from the ideal response of a perfect delay by 2*10' 5 . 
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Daubechies' 14-point generating filter deviates from an perfect delay by 1.4*10’ n . This 
represents six orders of magnitude in performance improvement over the time-domain 
performance of the QMF bank considered in Section III.C. 



Figure 5.7 indicates similar results for the frequency-domain performance of a QMF 
bank constructed from Daubechies' 14-point generating filter. The desirability of a 
power-complementary pair of analysis and synthesis filters was addressed in Chapter III. The 

QMF bank of Figure 3.15 deviates from being a perfectly power complementary pair by 10' 3 dB. 
A QMF bank constructed from Daubechies' 14-point generating function deviates from being a 
perfect power complementary pair by 1.2* 10' 1J dB. This represents eight orders of magnitude in 
improvement. Consequently, Daubechies’ generating function for wavelets and scaling functions 
supported on [0, 14] provide for outstanding QMF bank performance conforming very closely to 
the power complimentary and perfect reconstruction properties. 




CD 
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Figure 5. 7- Demonstration of the degree to which generating filter for Daubechies' orthonormal wavelet and 
scaling function on [0, 13] produces a power complementary pair for QMF bank purposes. The left-hand axis 
indicates the frequency responses for the filter and its modulated complement. The right-hand axis indicates the 
degree to which the filter pair forms a power complementary set. 
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D. 



DESIGN OF FILTERS FOR OMF BANKS 



[n the final Section of this chapter, procedures for the design of filters for filter bank 
applications will be examined. The steps for the design of a filter set for a two-channel QMF 
bank will be outlined. After that, the process of a filter design for a four-channel pseudo-QMF 
bank— a filter bank which only approximately satisfies the perfect reconstruction and power 
complementary properties— will be outlined. Procedures for the design of perfect reconstruction 
banks of an arbitrary number of channels exist [11, 34], but will not be addressed in the present 
study. 

To construct a two-channel filter bank with real-coefficient filters which possesses the 
perfect reconstruction property, it is necessary, as was shown by (3.20a), that the Z-transforms 
of the synthesis and analysis filters, both of order N, for a given channel k be related by [13] 

F k (z) = z' N H k ‘(l/z‘). (5.22) 

Equation (5.22) is equivalent to the notion that impulse response of the synthesis filter is equal 
to the time reversal of the impulse response for the analysis filter. Additionally, from (3.19), the 
equivalent transfer response T(z) of the filter bank structure must, to satisfy the perfect 
reconstruction criterion, be of the form 

T(z) =f [F 0 (z)H 0 (z) + F,(z)-H,(z)] = z' s (5.23) 

For convenience, an equivalent transfer response G k (z) for the branch of the filter bank 
corresponding to each channel can be defined: 

G k (z) = y-H k (z)-F k (z). (5.24) 
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Now, assume that G k (z) represents the transfer function for a half-band filter. If G k (z) 
satisfies the half-band filter property in the sense described by Mintzer [32], then the impulse 
response g k (n) corresponding to G k (z) satisfies [13] 

(g k ) i2 (n) = 5 Nn (5.25) 

where N represents the orders of the filters h k (n) and f k (n). Therefore, by (5.25), the only 

non-zero, even-numbered element of the impulse response of g k (n) is the sample g k (N). 
Continuing, the objective is to satisfy (5.23), or equivalently, 

T(z) = G 0 (z) + Gj(z) = z' N . (5.26) 

One way to satisfy (5.26) is to specify that the odd-numbered elements of g ( (n) each be of the 

opposite sign, or 

g t (2-n+ l) = : g 0 (2-n+ 1). (5.27a) 

Because of (5.25), (5.27a) is equivalent to 

gj(n) = (- 1 ) n ‘g 0 ( n ). (5.27b) 

But, (5.27b) is equivalent to 

G,(z) = G 0 (-z). (5.27c) 

From the preceding development, it may be inferred that H k (z) and F k (z) are 
equal-ordered spectral factors of a halfband filter G k (z). Furthermore, g,(n) is simply a 
modulated version of g^n). Therefore, the frequency response G^e* “) consists simply of the 
frequency response of G 0 (e l “) shifted in the frequency domain by n. Finally, since G k (z) is a 
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polynomial with two equal-ordered spectral factors which are simply the reversals of each other, 
certain constraints on the locations of the zeros of G k (z) arise. First, for each zero z q of G k (z) on 
the real axis corresponding to H k (z), there must be a reciprocal zero 1/z on the real axis 
corresponding to F k (z). Secondly, since it has been specified that the filter impulse responses 
h k (n) and f k (n) be strictly real, the zeros of G k (z) must occur in complex-conjugate pairs. 
Furthermore, for every zero z q located inside the unit circle corresponding to H k (z), there must 

be a corresponding reciprocal zero 1/z corresponding to F k (z) located outside the unit circle. 
Consequently, all zeros not on the unit circle or the real axis must occur in complex conjugate 
quadruples. Finally, in order for H k (z) and F k (z) to be of equal order, for each 
complex-conjugate pair of zeros z p and z q located on the unit circle corresponding to H k (z), an 
identical set corresponding to F k (z) must occur in the exact same location. In other words, the 
zeros on the unit circle for G k (z) must occur as double zeros. Consequently, G k (z) must be a 
filter with a strictly real, non-zero, symmetric, zero-phase filter. 

The requirements for the locations of the zeros for G k (z) suggest a design technique for 
filters for perfect reconstruction filter banks. The first step is to design a zero-phase, symmetric 
half-band, lowpass filter G 0 (z). This may be accomplished by any standard method such as the 
McClellan-Parks routine or by the window method [12]. Such a filter will be symmetric and 
will satisfy the half-band filter criterion (5.25). In order for G 0 (z) to be strictly non-zero, the 
maximum stopband ripple e R must next be measured. If this ripple e R is added to the central 
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element of the filter impulse response g 0 (n), the result will be an impulse response of a modified 
filter 



g 0R (n) = g 0 ( n ) + %-5( n - N) (5.28a) 

with a corresponding transfer function 

= G 0 (z) + (5.28b) 

Because of the addition of the term e R the frequency response of G 0R (z) will be strictly non-zero. 




Figure 5.<$-lmpulse response of 30-point, lowpass filter designed for QMF bank. 



The next step requires the rooting of G 0R (z). The zeros of G 0R (z) must be grouped into 
three subsets: zeros located inside of the unit circle, zeros located on the unit circle, and zeros 
located outside of the unit circle. The zeros inside of the unit circle belong strictly to H 0 (z) and 

the zeros outside of the unit circle belong strictly to F 0 (z). The zeros on the unit circle should all 
be double zeros occurring in sets of complex conjugate pairs. One set of each double complex 
conjugate pairs belongs to H 0 (z) while the other belongs to F 0 (z). If the zeros corresponding to 
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H 0 (z) are re-expanded into a polynomial, the coefficients of the result should be strictly real and 



belong to h 0 (n). 



As an example, a 30-point QMF filter was designed whose impulse response is plotted in 
Figure 5.8. The halfband filter was generated from applying a Kaiser window to a 59-point sine 
function 

sin I -S- (n-N)) 

g°(n) = — 



After adding the peak stopband error e R to impulse response element go(30), the corresponding 
polynomial was rooted and zeros grouped as was described above, resulting in a polynomial 
whose coefficients represent the impulse response h 0 (n) of a QMF-bank lowpass analysis filter. 




Sample sequence number, n 



C 

to 



Figure 5.9--Time-domain demonstration of performance of QMF bank constructed from the 30-point lowpass 
filter designed by spectral factorization of a half-band filter produced by the window method. Left-hand axis 
indicates equivalent impulse response while right-hand axis indicates deviation t(n) - 5(n - N Max ) from ideal 

perfect reconstruction performance. 



148 




Figure 5 . / (^-Demonstration of the degree to which 30-point lowpass filter designed by spectral factorization of a 
half-band filter produced by the window method produces a power complementary pair for QMF bank purposes. 
The left-hand axis indicates the frequency responses for the filter and its modulated complement. The right-hand 
axis indicates the degree to which the filter pair deviates from a power complementary set. 

Figures 5.9 and 5.10 indicate the results of application to the filter the analysis methods 
demonstrated in Chapter III for filter banks. Figure 5.9 plots the equivalent transfer response of 

a filter bank constructed from h 0 (n), and its deviation from the ideal response of an ideal delay. 
As can be seen from Figure 5.9, this method works quite effectively. The peak deviation from 
an ideal delay was less than 3x 10' 10 . This result compares favorably with the response for the 
generating filter for Daubechies' wavelets examined earlier. The power complementary 
properties are demonstrated in Figure 5.10. This filter did not quite perform as well as the 
Daubechies’ wavelet generating filter performed. The peak deviation from a response for a 

power complementary pair was approximately 6* 10' 9 compared with 1.2* 10' 11 for the 
Daubechies filter. This result is also slightly inferior to the performance of Daubechies' wavelet 
generating function by roughly two orders of magnitude, but still significantly superior to the 
QMF bank filters considered in Chapter III. 
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Finally, within the context of a study addressing wavelet transforms the question arises as 
to the use of perfect reconstruction QMF filters for generating wavelets. First of all, the filter 

impulse response h 0 (n) does not strictly satisfy (5.1 1). Specifically, 

Z h 2 (k) = -f + 7.5 x 10- 6 . 

k 

Consequently, any scaling functions or wavelets generated by h 0 (n) would not be orthonormal. 
Additionally, if a two-scale difference equation for h 0 (n) is set up in the form of (5.18), the 

nearest eigenvalue to unity for the double-right-shift matrix of coefficients of h 0 (n) is 1.0052. 

As a result, only an approximate system of a two-scale difference equation can be constructed 
from h 0 (n). Based on the preceding two measures it may be concluded that although, as was 




Figure 5. 1 /-Plots of "pseudo-scaling function" and "pseudo-wavelet" generated from four-stage dyadic expansion 
of 30-point filter designed from spectral factorization of a half-band filter. 

shown in Chapter IV, wavelet decompositions are exactly equivalent to structures of cascaded, 
perfect reconstruction QMF banks, the converse is not true. In fact, in the present example, the 
perfect reconstruction filter bank only approximately satisfies the properties described for 
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orthonormal wavelet decompositions. Nevertheless, if these issues are ignored and a dyadic 
expansion of h 0 (n) is performed, the resultant functions, after four iterations, appear as plotted in 
5.1 1. The functions possess obvious smoothness and present appearances similar to those of 
Daubechies' wavelets. 

To generate a filter bank for a four-subchannel multiresolution decomposition, a 
pseudo-QMF bank was designed using a spectral factorization technique similar to that used for 
the preceding example. The strategy was to approximate the power complementary property. 

Using the McClellan-Parks algorithm, two filters h 0 (n) and h,(n) were designed whose 
passbands covered the regions [0, jt/4] and [jt/4, tt/2], respectively. The upper half of the 
frequency spectrum below the Nyquist frequency was covered by modulating the first two 
filters: 

h 2 (n) = (-l) n -h,(n) 
and 

h 3 (n) = (-l) n h 0 (n). 

The design of the first two filters involved an iterative process in which, at each stage, the 
closeness of the filter set to the power complementary process was checked. At each step, the 
sum of the frequency responses was adjusted by modifying the boundaries to the passband and 
stopband for each filter. Tables 5.1a-c document the iterations. Each column of Tables 5.1a and 
5.1b indicates the frequency at which the desired gain in the top row is to be met for the filter to 
be factorized spectrally. Table 5.1c indicates, at each iteration, the peak deviation in the 
transition bands between each filter. Based on the deviation in the transition bands between 
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h 0 (n) and h,(n) and between h 2 (n) and h 3 (n), the filter indicated by h^n) was not modified after 
the third iteration. Eight iterations were necessary to obtain acceptable results for the transition 



band between h,(n) and h,(n). 

TABLE 5. L4--RECORD OF ITERATIVE ADJUSTMENTS TO PASSBAND AND STOPBAND FREQUENCY 



Target Gain 


1 


1 


10’ 8 


icr 8 


Frequency 


Iteration 1 


0 


tt/8 


5-71/16 


7t 


Iteration 2 


0 


3-71/16 


5-71/16 


K 


Iteration 3 


0 


21-71/128 


5-7T/16 


K 



TABLE 5 . /^-RECORD OF ITERATIVE ADJUSTMENTS TO PASSBAND AND STOPBAND FREQUENCY 



Target Gain 


10' 8 


10' 8 


1 


1 


10' 8 


10‘ 8 


Frequency 


Iteration 1 


0 


3-71/16 


5-7t/16 


7-71/16 


9-7T/16 


K 


Iteration 2 


0 


5-7T/32 


5-7t/16 


7-71/16 


19-7t/32 


X 


Iteration 3 


0 


19-7T/128 


2l-7t/64 


27-TT/64 


19-7t/32 


X 


Iteration 4 


0 


39-71/256 


41-71/128 


26-7t/64 


37-7t/64 


X 


Iteration 5 


0 


39-71/256 


83-7T/256 


53-7t/l28 


75-7T/128 


X 


Iteration 6 


0 


39-71/256 


83-7t/256 


105-7T/256 


151-71/256 


K 


Iteration 7 


0 


39-7t/256 


83-71/256 


52-7T/128 


151-71/256 


K 


Iteration 8 


0 


39-7C/256 


83-71/256 


103-7t/256 


15 1 -7C/256 


K 



Upon review of Table 5.1c and of Figures 5.13 and 5.14 it becomes obvious that the 
performance obtained from the four-channel filter bank is nowhere near what was obtained by 
any filterbank considered at any other time during this study. The peak delay of the equivalent 

impulse response for the filter bank was of the order of magnitude 5 X 10°. Furthermore, the 
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peak deviation from the power complementary standard could only be reduced to tenths of 
decibels. 

TABLE 5. 1 C-TRANSITION BAND DEVIATION FROM POWER COMPLEMENTARY, BY ITERATION, 



FOR DESIGN OF FOUR-CHANNEL PSEUDO-QMF FILTER BANK. 





Peak deviation in transition 

bands between h 0 (n) and h,(n) Peak deviation in transition 
and between h,(n) and h,(n) band between h^n) and h 2 (n) 


Iteration 1 
Iteration 2 
Iteration 3 
Iteration 4 
Iteration 5 
Iteration 6 
Iteration 7 
Iteration 8 


Excessive 


Excessive 


0.22 


0.38 


0.07 


0.23 


0.08 


0.2 


0.02 


0.2 


0.02 


0.08 


0.02 ' 


0.06 


0.02 


0.03 




Figure 5. 12 - Impulse responses h 0 (n) and h,(n) of filters designed by spectral factorization of fourth-band filters 

designed using McClellan-Parks algorithm. 
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Figure 5. 1 J-Time-domain demonstration of performance of four-channel QMF bank constructed from filters 
designed from spectral factorization of fourth-band McClellan-Parks filters. Left-hand axis indicates equivalent 
impulse response while right-hand axis indicates deviation t(n) - S(n - N Max ) from ideal perfect reconstruction 

performance. 
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Figure 5. 14 - Demonstration of the degree to which four-channel, pseudo-QMF bank designed by spectral 
factorization of a McClellan-Parks fourth-band filters produces a power complementary pair for QMF bank 
purposes. The left-hand axis indicates the frequency responses for the filter and its modulated complement. The 
right-hand axis indicates the degree to which the filter pair deviates from a power complementary set. 



154 



VI. EVALUATION OF MULTIRESOLUTION TECHNIQUES FOR 

DETECTION APPLICATIONS 

A. INTRODUCTION 

In this chapter, additional performance characteristics of the multiresolution structures 
developed in Chapter IV will be examined within the context of evaluation of suitability for 
detection applications. In Section B, the computational efficiency of multiresolution structures 
will be considered. In the Section C, receiver operating characteristics will be plotted for 
various multiresolution structures applied to the test sequences generated by (4.74b) and (4.74a). 
Finally, Section D will demonstrate resolution of scale within the context of steady-state 
harmonics and chirped signals. 

The basis for evaluation of multiresolution performance will be the spectrogram, one of 
the most commonly used techniques for detection and spectral estimation. Figures 6. la and 
6.1b present the spectrogram decomposition of the two 256-point test sequences generated by 
(4.74a) and (4.74b), respectively. The plots were generated by plotting the frequency spectrum 
below the Nyquist frequency for a sequence windowed with a 128-point Hanning window 
shifted in increments of two samples. This representation was selected because of its frequency 
resolution and because it presents a similar number of points compared to the multiresolution 
decompositions presented in Chapter IV. 

In Figure 6.1a, the lack of frequency resolution in the extreme lower range of the 
frequency spectrum is apparent. For lowpass test sequence (4.74a), the spectral components at 
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digital frequencies 2 -jt 2/512 and 2-71-9/512 are not clearly resolved. The component at 
2-7C-2/5 1 2, which appears at spectral bin location 0.5, because of its relatively greater power, 
spreads into adjacent spectral bins. Consequently the component at 2-71-9/512, which coincides 
with spectrogram bin 2.25, is partially masked by leakage from the component at 2 -Jt- 2/5 12. The 
only indication of the component at 2-71-9/512 is evident on the leading edge of the surface plot. 




Figure <5./a--Spectogram representation of lowpass test sequence (4.74a). The plot was generated using a 
128-point Hanning window shifted in increments of two samples. Only the half of the frequency spectrum below 

the Nyquist frequency is plotted. 

On the other hand, Figure 6.1b provides an obvious indication of the advantages of 
constant spectral bin density: The three harmonics of (4.74a) are quite easily resolved. The 
plots have the additional advantage of a smoother presentation. Spectrogram decompositions 
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smooth out the phase information while multiresolution presentations indicate individual 
oscillation cycles of the sequences decomposed. 
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Figure 6. 76--Spectogram representation of highpass test sequence (4.74b). Tlie plot was generated using a 
128-point Hanning window shifted in increments of two samples. Only the half of the frequency spectnim below 

the Nyquist frequency is plotted. 

B. COMPUTATIONAL EFFICIENCY OF SIGNAL ANALYSIS TECHNIQUES 

Computational efficiency comprises one of the principal advantages ascribed to Fast 
Fourier Transform (FFT) methods of digital signal processing. The required number of 
multiplications provides a common number for evaluating digital processing techniques. In the 
case of FFT decomposition of a sequence of length M, the required number of complex 

multiplications N x is expressed by 

N x = flog 2 (M). (6.1) 
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Furthermore, application of the window sequence to each segment of a sequence of length L 
requires M multiplications. Finally, if the window is shifted in increments of k, then the FFT 
will be applied to -^+1 sequences. Therefore, the total number of multiplications necessary 
to generate a windowed FFT decomposition is 

N * = (^ +1 )T- log 2 (M) - (6.2a) 

In general, the multiplications performed in FFT decomposition represent complex 
multiplications. Evaluation of a complex multiplication requires four real multiplications and 
two additions. Consequently, assuming that real data and window sequences are used, the 
number of multiplications for the FFT decomposition is actually 

N*= 2^+1 )m 2 log 2 (M). (6:2b) 

For each of the decompositions plotted in Figures 6.1a and 6.1b, the sequence length L = 256, 
and the window length M = 128. Furthermore, the window was shifted in increments of k = 2 
samples. Consequently, using (6.2b), the number of real multiplications N x required to perform 
the decomposition was N x = 14,909,0440. 

For the multiresolution decompositions of Chapter IV, Section E, FFT-type radix-2 
algorithms are not available. Consequently, each stage of decomposition requires the 

convolution of two real sequences. Therefore, for decomposition of sequences of length L 0 , 

using filters of length M c for decomposition into two primary channels, the required number of 
real multiplications for the first stage is 

N< = 2-M c -L 0 . 
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After convolution by the analysis filters of each stage, the approximation and detail sequences 
are both decimated, providing sequences of length 

» Li)+Mc~l 

Ll = 2 ' 

By induction, the length L p of the approximation and detail sequences at stage p is given by the 
geometric series 

L p = t£ + (M c - 1) • Z ± 

1=1 

= -^r + (M c - 1) ■ (1 -j?) ■ (6.3) 

= tj- • ^L 0 -(M c - 1) j + (M C - 1) 

Since, at each stage, each of two sequences is convolved with a filter of length M c , the total 

number of multiplications N c necessary to perform a K-stage multiresolution decomposition into 
primary channels 

Nc=2-M c - z‘l p . (6.4) 

P =o 

Substitution of (6.3) into (6.4) produces 

Nc=2-K-(Mc-l) + 2-M c - (Lo-(Mc-l)j- z'^. (6.5) 

Finally, evaluating the geometric series in (6.5) produces 

N c = 2- K- (M c - 1) + 4- M c • ( L 0 -(Me - 1) ) • (1 --±r). (6.6) 

Equation (6.6) covers all of the multiplications necessary to perform the decomposition into 
primary channels in accordance with block diagram 4.19. 
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If the primary detail channels are each to be decomposed into J subchannels using a 
subchannel filter bank constructed from filter of length M s , additional computation is necessary. 
First, if a sequence of length L p is convolved with a filter of length M s and subsequently 
decimated by a factor of J, the length L r (p) of the resultant sequence will be 

L r ,p, = }(L P + M S -1). 

The number of multiplications N r (p) to perform this signal reduction at stage p is expressed by 

N r ' p) = M S -L P , 

which by substitution of (6.3) becomes 

N[ p) = Ms • • (l 0 - (M c - 1 )) + (M c - 1 )) . 

If there are a total of J subchannels per stage and a total of K stages, then the total number of 
multiplications N r required to perform all subchannel decomposition operations is evaluated as 

N r = J- I Nr P) 

p = l 

= J.K-M,-(M c -l) + J-M I ^L 0 -(M c -l)j- (6.7) 

= J- K- Ms • (M c - 1) + J-M S • (Lo-(M c - 1)) • (1 -■£) 

If the decomposed sequence of length L r lp) is expanded by a factor of J and subsequently 
convolved with a filter of length M c , the length L e <p) of the resultant sequence is expressed by 

L e (p) = J-L r ,p) + M s - 1 = L p + 2-(M s - 1 ). 

The number of multiplications N e <p) necessary to perform this operation for J subchannels at 
decomposition stage p is therefore 
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N, 11 ” = JM, L,"” = J M, (L„ + 2 (M, - 0). 

Furthermore, the total number of multiplications necessary to expand all subchannel 
decompositions is expressed as 



N e = J • M s • I L« p) = J • M s • X !’^-fL, J -(Mc-l)) + (M c -l)j. 



K L7' = J • Ms • X | 5? 

p=i p=iL 2 

Substituting (6.3) into (6.8) and evaluating the series produces 

K 



( 6 . 8 ) 



= J-M s - Z 

P = lL 



^■•(Lo-(Mc-l) I +(M C - l) + 2 • (M s - 1) 



1 
J 

= J • K • M s • |^(M C - 1) + 2 • (M s - 1) J + J • M s • |^Lo -(M c - 1) j • I j? 

= J • K • Ms • |(M C - 1) + 2 • (Ms - 1)] + J • Ms • (l 0 - (M c - 1) j • ( 1 - ^) 



(6.9) 



To calculate the number of multiplications necessary to re-expand a detail component 
sequence reconstructed from decomposition into subchannels, the notation L q lp) will be 
employed to indicate the length of the sequence at a stage of interest. The superscript (p) 
indicates the number of the stage of decomposition at which the detail component sequence was 
extracted. The subscript denotes the stage of reconstruction. The stages of reconstruction are 

labeled in decreasing order from p-1 to 0 for a component extracted at the p 111 stage of 
decomposition. At stage p, the length of the detail sequence, after reconstruction from 
subchannel decomposition, is given by L e (p) = L p (p) . Now, upon expanding the detail sequence of 
length L p lp) and convolving it with the channel synthesis filter of length M c , at stage p-1, the 
expanded sequence which is obtained has length 



L , (p) = 2-L (p) + M c - 1 = 2 -L + 4-(M s - 1) + M c - 1. 



( 6 . 10 ) 
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Inductively, after q stages of decimation followed by convolution with the appropriate channel's 
synthesis filters, (6.10) becomes 

l P T = 2 q -(L p + 2-(M s -l)) + (M c -U- zV 

= 2 q - (L p + 2-(M s -l)j + (M c -l)-(2 q -l) . (6.11) 

= 2 q • |^Lp + 2 • (M s - 1) + (M c - 1) j - (M c - 1) 

Now, substituting L p from (6.3) into (6. 1 1) produces 

Lp_q = 2 q -^-^Lo-(Mc-l)) + (Mc-l) + 2-(Ms-l) + (M c -l)j-(M c -l) . 

(6.12) 

Consequently, the total number of multiplications required for a full expansion by channel of a 
decomposed sequence can be evaluated by summing up the contributions from each convolution 

performed in each channel. It stands to reason, that the number of multiplications 

involved in expansion of J detail component sequences (one detail component sequence for each 

of J subchannels) of length L p lp> , at stage p, to obtain a sequence of length L p ., lpl ,at stage p - 1, is 
indicated by 

Inductively, therefore, the total number of multiplications N 0 (p) necessary to completely 
expand the J decomposed sequences in channel p is obtained from 

Nj p) = J-M c - I Lp-i P) . (6.13) 
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Therefore, to obtain N 0 (p, ( it is necessary to substitute (6.12) into (6.13) and evaluate the series, 
which produces 



Nf> 



= J Me 



-T- Lo-(M c -l) +2 • (M c - 1) + 2 ■ (M s - 1) 



I (2 q ) -p ■ J ■ M c • (M« - 1) 

q=l 



= 2 • J • M c 



1 V (Lo-(M c -l) ]+2 (M c -l) + 2 (M s -1) 



(2 P - 1) - p ■ J • M c (M e - 1) 



= 2 • J • Me ■ ( Lo - 3 ■ (Me - 1) - 2 • (M, - 1 ) J + 2 • J • M c -[2 ■ (M« - 1) + 2 • (M, - 1 ) J • 2 P 

-2-J-Me- (Lo-(Me-l)j--^-p- J Me-(Me-l) 

(6.14) 



And lastly, to obtain N 0 , the total number of multiplications to compute the full 

expansion by channel, it is necessary to sum the contributions N 0 (p) from each channel. This 
sum appears as 



No = 2 • J • K • M c • ^L 0 -3 • (M c - 1) - 2 • (M s - 1) j 

+ 2 • J - Me • ^2 • (M c - 1) + 2- (M s - 1) j • I 2? 
-2 - J-M c -(l 0 -(M c -1)) - £ £ 

-J-Me (Mc-l)- £ p 



(6.15) 



The first series in (6.15) are simply the geometric series evaluated previously in this 
development. The last series can be solved using Bernoulli polynomials [40] which produce 

j, p =H KI+ 4 

Substituting each of the evaluated series into (6.15) and collecting similar terms produces, for 

N 0 , 



163 



No = 2 • K • J • M c • ^L 0 -3 • (M c - 1) - 2 • (M s - 1) j 

+ 4-J-M c -(V(Me-l) + 2-(M s -l) 

-2-J-M c -(Lo-(M c -1))-(1-^) 

-Y-J-Me-(Me-l)-^K 2 +Kj 




Equation (6. 16) applies whether or not subchannel decomposition is performed. If 
subchannel is not performed, then the detail sequence at each stage is decomposed into J = 1 

subchannels and the filter g 0 (n) used for subchannel decomposition possesses as its impulse 



response g 0 (n) = 5 0 n where 5 0 n is Kronecker's delta function. Consequently, because of how 

g 0 (n) is defined, M s =l and, in (6.16), each of the terms in involving M s - 1 become zero. 

TABLE 6. 7-COMPARISON OF NUMBER OF MULTIPLICATIONS NECESSARY FOR MULTIRESOLUTION 



AND SPECTROGRAM SIGNAL REPRESENTATIONS. 





Primary Primary 

Channel Channel Subchannel 

Decomposition Re-expansion - Decomposition 


Subchannel 

Re-expansion 


Total 


Two-Channel, 

Zero-Subchannel 


13,762 


406,558 


0 


0 


420,321 


Two-Channel, 

1 Two-Subchannel 


13,762 


1,544,028 


9,689 


15,513 


1,582,902 


Two-Channel, 

Three-Subchannel 


13,762 


2,400,378 


15,572 


25,625 


2,455,274 


Two-Channel, 

Four-Subchannel 


13,762 


4,887,224 


41,526 


97,206 


5,039,634 


Spectrogram 










14,909,440 











Summarizing, there are four contributors which must be evaluated to calculate the 
number of multiplications N 0 necessary for a full channel by expansion of a decomposed 
sequence. The basic decomposition, N c , is provided by (6.6). The decomposition into 
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subchannels, N r is calculated using (6.7). To re-expand subchannel decompositions, the number 
of multiplications is indicated by (6.9). And finally, for the full expansion by channel, the 
number of calculations is provided by (6. 15). For some of the multiresolution structures 
demonstrated in Section IV.E, the number of multiplications is tabulated in Table 6.1. 

In considering Table 6.1, it must be noted that the decompositions listed do not represent 
the optimum representation. In the case of the multiresolution decompositions, the essential 
information about the decomposed signal is contained in the primary channel and subchannel. 
The re-expansion by channel merely provides a manner of displaying the information which is 
more convenient than plotting coefficients in a non-uniform lattice of the type of Figure 4.8. 

In the case of the spectrogram, shifting a 128-point window by increments of two 
samples provides what is likely excessive overlap. This degree of overlap was used in an 
attempt to capture, to the greatest extent possible, the time-varying features of the test sequences 
on which decomposition was performed. If, for instance, the 128-point window was shifted in 

increments of 64 samples to provide for a 50% overlap, the number of multiplications N x would 

be N y =688,128. This overlap, however, would provide poorer indication of the non-stationary 
features of the sequences for which decomposition was performed. A 50% overlap with a 
window length of 128 used to decompose a sequence of length 256 would provide for three time 
shifts. 

With regard to the multiresolution decompositions, excluding the calculation used to 
re-expand the components of decomposition, a significant advantage in computational efficiency 
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is evident. Because of the factor-of-two decimation performed at each stage of decomposition, 
the sequence length at each stage is one half the length of the sequence of the previous stage. 
Consequently, the number of calculations necessary to perform the calculation is related to a 
geometric series in one half. 

C. PERFORMANCE OF DECOMPOSITIONS Pi THE PRESENCE OF NOISE 

Many signal processing techniques are employed with the intent of extracting a signal 

from noise in which it is embedded. The spectrogram, one of the most commonly used signal 
processing tools, represents an attempt to localize spectral components of a signal. If the noise 
in which the signal is imbedded can be characterized as "white," possessing a flat power spectral 
density, a signal which can be localized in a spectral sense will be detectable in greater noise 
intensity than a signal which is not spectrally localized. Consequently, when attempting to 
detect a signal by means of spectral localization, the effective signal-noise ratio can be reduced 
to the ratio of signal power to the power of only the noise contained within the spectral region of 
interest. 

To detect a time-varying process, the detection process becomes a two-dimensional 
problem, requiring localization in both frequency and time. Within the context of 
time-frequency localization, multiresolution techniques present the potential for an alternative 
method of detection. Many of the advantages of multiresolution ttansforms which lend those 
decompositions to signal decomposition also provide the basis for extraction of embedded 
signals from noise. Multiresolution transforms provide good time resolution at high frequencies 
and frequency resolution at lower frequencies. Furthermore, if it is sought to detect a class of 
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processes for which bandwidth is proportional to the center frequency of the process, 
multiresolution techniques may provide the optimum detection technique. 

If a sequence s(n) is embedded in additive, white, Gaussian noise q(n), the resultant 
sequence 

x(n) = s(n) + q(n) (6.17) 

will also be a Gaussian random variable. If the noise process has rj(n) a mean of zero then the 
mean of the signal plus noise x(n) is indicated by 

x(n) = s(n). (6.18) 

Furthermore, the variance of x(n) equals the variance of q(n): 

c x 2 = a n 2 . (6.19) 

In Chapter IV it was demonstrated that decomposition of a sequence by multiresolution 
structures constructed from cascaded perfect reconstruction quadrature mirror filter banks is 
equivalent to projection of the sequence to be decomposed upon a space of orthonormal vectors. 
Decomposition by perfect reconstruction Filter banks into subchannels, furthermore, can be 

interpreted as subdivision of the detail vector space W m , using the notation from Chapter IV, 

into another set of complementary set of subspaces { W m ' k) } k . If, however, the filter banks do 
not satisfy the perfect reconstruction criteria, orthonormality of the equivalent projection 
operation is not guaranteed. 

The the wavelet coefficients {b m n } nsZ obtained from decomposing a sequence x(n) 
were defined in Section IV. C as 
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X, Vl/ m> n 



( 6 . 20 ) 



bm, n — 

Now, if x(n) is a Gaussian random variable, then b m n will also be a Gaussian random variable. 

Furthermore, if the set of basis functions {vj f m n } is orthonormal, then the wavelet coefficients 

b m n will be linearly independent. To obtain the expected value of b m n , substitution of (6.17) 
into (6.20) produces 

b m,n = ( S >V m ,„) + (^¥ m J. (6.21) 

Assuming the T| has a mean of zero, then the mean of the wavelet coefficient b m n for a signal 
s(n) embedded in noise is simply the value of the wavelet coefficient evaluated in the absence of 
noise. 

To obtain the variance of b m n , the squared mean is first 

£{b m .„} = N, M/ m . J 2 + C{(n, (6.22) 

The expectation on the right-hand side of (6.22) is evaluated as 

S{(T¥ m ,n) 2 } - £{JJ fi(u)ri(v)v m , n (u)Vm.n(v) dudvj. (6.23) 

The expectation operator on the right-hand side of (6.23) can be passed through the double 
integral to operate on the terms r|(u)-r|(v) producing 

£{fi(u)-n(v)} =a n 2 -5(u- v) (6.24) 

where 5(t) denotes the Dirac delta function. By the integral sifting property of the Dirac delta 
function, (6.24) reduces to 

^Ol’VmJ 2 } = V(¥ m .„.¥ m ,n) = V (625) 
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The wavelet coefficients, therefore, simply comprise Gaussian random variables whose mean is 
(s, \j/ m n ) with a variance of o^ 2 . 

If detection of a known sequence in the presence of noise is required, one way to generate 
a decision statistic is to sum the squared magnitudes of all of the nodes in the transform space 
whose powers exceed a specified threshold. In performing multiresolution decomposition of a 
sequence, the essential information about a signal is contained in a linearly independent array 

Dj k of non-uniform density such as that in figure 4.8. In detecting a known sequence in the 
presence of noise, it is conveniently possible to consider only the decomposition nodes in which 
a significant portion of the signal energy is represented. Mathematically, given a multiresolution 

decomposition represented by a lattice of points Dj k , the decision statistic can, equivalently, be 
constructed which excludes nodes contained in k which are negligible in value. The subset of 
D ; k of which the decision statistic is comprised can be limited to those nodes whose addresses 
are indicated by 1C 



O’, k) : 


l l 2 


| 2 


D,k >p-max 


Dj.k 1 




1 J> 1 j.k 



Figures 6.2a, 6.2b, 6.3a and 6.3b present the results for statistics constructed from all 
nodes defined by (6.26) where p = 1/100. Table 6.2 indicates the mean square reconstruction 
error (3.32) which results from the value selected for p. The information in Table 6.2 is 
intended as an indication of the consequences incurred from considering only the nodes in D k 
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where k is defined by (6.26). The decision statistics for the periodogram were generated in the 
same manner as for the multiresolution decompositions. 



TABLE <5.2— NORMALIZED MEAN-SQUARE RECONSTRUCTION ERROR RESULTING FROM 
TRUNCATED DECOMPOSITION. RECONSTRUCTION PERFORMED FROM SUBLATTICE DEFINED BY 
(6.26) WITH p = l/l 00. EIGHT DECOMPOSITION STAGES WERE PERFORMED IN THE TWO-CHANNEL 
CASES WHILE FIVE-STAGE DECOMPOSITIONS ARE INDICATED FOR THE THREE-CHANNEL CASES. 





Test sequence (4.74a) 


Test Sequence (4.74b) 


Two-channel, 

zero-subchannel 


0.0065 


0.0121 


Two-channel, 

two-subchannel 


0.0121 


0.0244 1 


Two-channei, 

three-subchannel 


0.0046 


0.0209 


Two-channel, 

four-subchannel 


0.0288 


0.1578 


Three-channel, 

zero-subchannel 


0.01 


0.01 




Figure 6.2a — Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 
low-frequency transient (4.74a) embedded in noise with a -3 dB signal-noise ratio. Characteristics generated from 

500-realizations of 256-point random noise vector. 

For each technique indicated, 500 realizations of a 256-point, Gaussian random vector 



were generated. For each multiresolution technique reflected in Figures 6.2a and 6.2b, prior 
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generating the receiver operating characteristics, decision statistics were generated for separate 
noise realizations alone. The receiver operator characteristics for noise alone consisted of 
unity-slope straight lines. The decompositions were then performed for the test sequences 
(4.74a) and (4.74b) with additive noise and then for the noise vector alone. The decision 
statistics for the receiver operator characteristics consisted of all nodes in the composition whose 
magnitudes were greater than the one-percent threshold described above. 
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Figure <5.26--Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 
high-frequency transient (4.74b) embedded in noise with a -3 dB signal-noise ratio. Characteristics generated 

from 500-realizations of 256-point random noise vector. 

The receiver operating characteristics were, therefore, constructed from central and 
non-central chi-square distributions. The sum of the squares of the selected nodes of the noise 
decomposition, since they were generated from zero-mean noise, constitutes a chi-square 
distribution. Summing the selected nodes of the decomposition of signal plus noise produces a 
non-central chi-square distribution. 



In both Figures 6.2a and 6.2b, the spectrogram provided greater robustness than the 



multiresolution techniques. With the exception of the three-channel, zero-subchannel 
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decomposition applied to highpass test sequence (4.74b) in -3 dB signal-noise ratio, the 
differences were not, however, dramatic. Moreover, for test sequence (4.74b), the 
multiresolution curves were, with the named exception, closer to the performance of the 
spectrogram than for lowpass test sequence (4.74a). The greater disparity in the case of highpass 
test sequence (4.74b) may be attributable to the difference in resolution between the 
multiresolution and spectrogram techniques in the high-frequency regions of the spectrum. 

It is interesting to observe that, in Figure 6.2a, the best results among the multiresolution 
decompositions were obtained from the two-channel, zero-subchannel and two-channel, 
two-subchannel algorithms. Those two algorithms were constructed strictly from 
perfect-reconstruction filter banks and, therefore, produced the lowest reconstruction error. 
Additionally, in the case of highpass test sequence (4.74b), the best performance was obtained 
from the two-channel, three-subchannel decomposition. The two-channel, four-subchannel 
decomposition provided, by a small margin, the worst of all results. 




Figure 6. 3a - Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 
low-frequency transient (4.74b) embedded in noise with a -6 dB signal-noise ratio. Characteristics generated from 

500-realizations of 256-point random noise vector. 
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Figure 6.3b— Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 
high-frequency transient (4.74a) embedded in noise with a -6 dB signal-noise ratio. Characteristics generated 

from 500-realizations of 256-point random noise vector. 

Results for a signal-noise ratio of -6 dB are presented in Figures 6.3a and 6.3b. In the 
case of the lower signal-noise ratio, all of the multiresolution methods outperformed the 
spectrogram by a small measure. In the case of lowpass test sequence (4.74a) the best results 
were provided, as with the -3 dB signal-noise ratio test, by the two-channel, two-subchannel and 
the two-channel, zero-subchannel, respectively. This trend may be related to the superior 
reconstruction error provided by these techniques. In the case of highpass test signal (4.74b), on 
the other hand, the best results came from the four-subchannel and three-subchannel 
decompositions. This may have occurred as a result of the resolution of the decompositions. 

For the low-frequency transient, (4.74a), the three-subchannel and four-subchannel may have 
over-resolved the signal. For the high-frequency transient, (4.74b), the two-subchannel and 
zero-subchannel decompositions may have under-resolved the signal. Consequently, in the first 
case, the signal energy, for the higher resolution decompositions, may have been distributed over 
a greater-than-optimum number of time-scale bins. In the second case, the lower resolution 



173 



decomposition may not have sufficiently localized the signal to isolate it from the effects of the 
noise. 

D. DECOMPOSITION OF STEADY - STATE HARMONICS AND FM CHIRPS 
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Figure d.4a--Time-frequeucy plot of evolution of constituents of first 256-point test signal constructed from one 
steady-state harmonic and a quadratically chirped component. 
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Figure 6. ^^-Time-frequency plot of evolution of constituents of second 256-point test signal constructed from 
one steady-state harmonic and a quadratically chirped component. 

The time-scale resolutions of multiresolution decomposition techniques were somewhat 



addressed in Section IV.E. However, in order to better understand the performance of these 



methods, a more controlled demonstration is desirable. In this section, results will be 
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demonstrated for a frequency-modulated, quadratic chirp converging to a steady-state harmonic 
component. The performance if the spectrogram method will first be illustrated. Then the 
results will be considered for two multiresolution techniques— the two-channel, four-subchannel 
and the two-channel, three-subchannel decompositions. 

Two test signals were generated for this demonstration. To illustrate the variations of 
resolution of scale provided by multiresolution decompositions in different regions of the 
frequency spectrum, it was desired to demonstrate the case in which the quadratic FM chirp 
converges from above to a steady-state harmonic in the lower half of the frequency spectrum and 
the case in which the FM chirp converges from below to a steady-state harmonic in the upper 
half of the frequency spectrum. Based on the experience of Chapter IV, it is expected that, for 
the first case, better resolution will be provided by the multiresolution techniques. 

The form for dynamic component of the test sequences is given by 



where the instantaneous frequency fft) possesses the form of a parabola whose vertex is, in the 
time-frequency plane, located at (t p f f ), the coordinates of the final time and frequency in the test 
region. The time-frequency parabola satisfying this criterion assumes the form of 



s(t) = cos(2 • rt • fit) dt) + cos(2 • Jt • f f • n) 



(6.27) 



f(t) = oc-(t - t f ) 2 + f f . 



(6.28) 



Evaluating the integral in the argument of the first cosine term on the right-hand side of (6.26) 



produces for the phase of the chirped component 




(6.29) 
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To uniquely define fit), it is necessary to specify the location of (t f , f f ), the veitex of the 
parabola, and the location of some other point on the curve in the time-frequency plane. If the 
time t 0 = 0 and frequency f 0 of the chirp at the beginning of the test sequence are specified, the 
coefficient a from (6.28) becomes 

Ct=^-r-. ( 6 . 3 ( 

! r 

Since 256-point test sequences are to be used, the final time t f = 256. It therefore remains to 
specify beginning and ending frequencies. 



For the first test sequence, a chirp which converges from below to a steady-state 
harmonic in the upper half of the frequency spectrum below the Nyquist frequency. The 
beginning and ending frequencies selected for the chirp were 

fo = 7k- f ‘ 

and • (6.31a) 

fr = -k'i 

For the second test sequence in which the quadratic chirp converged from above to a steady-state 
harmonic in the lower half of the frequency spectrum, the frequencies selected were 

fo = i-fs 

and • (6.31b) 

The ideal, time-frequency evolution of the two test signals are plotted in Figures 4.6a and 4.6b, 
respectively. 

The spectrogram decompositions of the two test signals are plotted in Figures 6.5a and 
6.5b for two purposes. The first reason is to verify that the desired affects were obtained in 
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construction of the test sequences. The second reason for the periodograms is to provide a 
standard with which the performance of the multiresolution decompositions can be compared. 
The spectrogram plots were generated in a manner identical to that used in Section A of the 
present chapter. The data sequences were windowed by a 128-point Hanning window which 
was shifted in two-sample increments. Only the half of the frequency spectrum below the 
Nyquist frequency was plotted. 




Figure 6. 5a--Spectrogram decomposition of first 256-point test sequence constructed from a steady-state 
harmonic component and a frequency-modulated quadratic chirp. The decomposition was obtained by windowing 
the data sequence with a 128-point Hanning window shifted in increments of two samples. Only the half of the 
frequency spectrum below the Nyquist frequency is plotted. 

Examining Figures 6.5a and 6.5b, the two components of each test sequence can be seen 
to converge as expected. The steady-state harmonic produces a much better resolved indication 
than the chirped component. This occurs because of the smearing of the chirp which results 
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from the frequency shift occurring during the integration time of each FFT. Both 
decompositions represent the data with uniform resolution throughout the frequency spectrum. 
Consequently, the time, for both sequences, at which the components cannot be resolved is the 

same for both sequences at approximately the 230 th sample. Consequently, from (6.28), the 



spectral resolution of the spectrogram decomposition is approximately 0.0043 f s . No attempt 
was made to optimize the window length with respect to resolution of the test signal 
components. The 128-point window provides good frequency resolution, however, because of 
the long integration time, the frequency of the chirped component has shifted significantly 
during each transform. 
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Figure 6. 5b - Spectrogram decomposition of second 256-point test sequence constructed from a steady-state 
harmonic component and a frequency-modulated quadratic chirp. The decomposition was obtained by windowing 
the data sequence with a 128-point Hanning window shifted in increments of two samples. Only the half of the 

frequency spectrum below the Nyquist frequency is plotted. 
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Figure 6. <5a--Two-channel, three-subchannel multiresolution decomposition of first 256-point test sequence 
constructed from steady-state harmonic component and quadratic chirp. Primary channel filters were based on 
Daubechies' orthonormal wavelet and scaling function on [0, 13]. Subchannel decomposition was performed with 

pseudo-QMF filter bank designed in [1 1], 

Examining Figures 6.5a and 6.5b, the two components of each test sequence can be seen 
to converge as expected. The steady-state harmonic produces a much better resolved indication 
than the chirped component. This occurs because of the smearing of the chirp which results 
from the frequency shift occurring during the integration time of each FFT. Both 
decompositions represent the data with uniform resolution throughout the frequency spectrum. 
Consequently, the time, for both sequences, at which the components cannot be resolved is the 
same for both sequences at approximately the 230 th sample. Consequently, from (6.28), the 
spectral resolution of the spectrogram decomposition is approximately 0.0043-f s . No attempt 
was made to optimize the window length with respect to resolution of the test signal 
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components. The 128-point window provides good frequency resolution, however, because of 
the long integration time, the frequency of the chirped component has shifted significantly 
during each transform. 
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Figure 6.6b— Two-channel, three-subchannel multiresolution decomposition of second 256-point test sequence 
constructed from steady-state harmonic component and quadratic chirp. Primary channel filters were based on 
Daubechies' orthonormal wavelet and scaling function on [0, 13]. Subchannel decomposition was performed with 

pseudo-QMF filter bank designed in [l 1]. 

The two-channel, three-subchannel decomposition of the second test sequence is plotted 
in Figure 6.6b. In the case of this plot, the chirp becomes unresolvable from the steady-state 

harmonic at approximately the 230 th sample. By (6.28), therefore, the frequency resolution is 

approximately 0.0043-f s . This number is quite close to the Figure for the spectrogram 
decomposition. Furthermore, the two test sequence components become unresolvable around a 
scale of four which coincides detail subchannels obtained from the fourth stage of 
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decomposition. If the final frequency f f had occurred at a lower frequency, the convergence 
point between the two components would have occurred in the detail subchannels of a later 
decomposition stage resulting in even greater resolution. 
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Figure 6. 7a-- Two-channel, four-subchannel multiresolution decomposition of first 256-point test sequence 
constructed from steady-state harmonic component and quadratic chirp. Primary channel filters were based on 
Daubechies’ orthonormal wavelet and scaling function on [0, 13]. Subchannel decomposition was performed with 

pseudo-QMF filter designed in Chapter V. 

Figures 6.6a and 6.6b present results of three-channel, four-subchannel multiresolution 
decomposition of the two test sequences. In Section IV.E, within the context of considering the 



decomposition of transients constructed from components of varying frequencies, the 
performance of the three-channel decomposition proved quite similar to that of the four-channel 
decomposition. This similarity appears in the case of decomposition of the present set of test 
sequences. In Figure 6.7a, the chirp and the harmonic become indistinguishable at 
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approximately sample number 108. This represents a slight improvement over the 
three-subchannel decomposition. In the case of Figure 6.7a, again using (6.28), the frequency 

resolution appears to be approximately 0.140f s . This represents a subtle improvement over the 
performance of the three-channel decomposition. Nevertheless, in this region of the frequency 
spectrum, the spectrogram continues to provide performance which is superior by a significant 
margin. 




Figure 6. 7b— Two-channel, four-subchannel multiresolution decomposition of second 256-point test sequence 
constructed from steady-state harmonic component and quadratic chirp. Primary channel filters were based on 
Daubechies' orthonormal wavelet and scaling function on [0, 13]. Subchannel decomposition was performed with 

pseudo-QMF filter designed in Chapter V. 

Finally, the performance decomposition of the second test sequence is plotted in Figure 
6.7b. The chirp and the harmonic converge at approximately the same location. This result is 
also consistent the observations of Section IV. E. When decomposing transient test signals, it 
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was observed that increasing from three to four subchannels does not dramatically improve 
resolution of scale at large scales. The primary from the increase is realized in the smaller 
scale— or higher frequency— regions. 
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VII. CONCLUSION 



A. SUMMARY OF RESULTS 

When compared with conventional periodogram techniques, the multiresolution techniques 
considered provided mixed results. In resolving signals at the extreme low end of the frequency 
spectrum below the Nyquist frequency, the multiresolution structures achieved their best results. 
In higher frequency regions, the multiresolution proved to be on par with or inferior to the 
periodogram. Overall, wavelet-like methods should be optimum for proportional bandwidth 
processes. The test signals employed performed marginally at demonstrating this strength in the 
upper regions of the frequency spectrum. 

With respect to other performance measures, the multiresolution techniques achieved, even when 
expanding fully by channel, decomposition of the test sequences with fewer calculations than the 
periodogram decomposition. In fairness to the latter method, the implementation demonstrated 
did not reflect any attempts to optimize the technique for the applications shown. However, in 
determining the number of calculations, the computational burden necessary to convert the 
resulting complex decomposition into a real representation was not considered either. 

The most interesting results were perhaps obtained from generation of the receiver operating 
characteristics. In the relatively high (-3 dB) signal-noise ratio demonstration the periodogram 
outperformed the multiresolution decomposition techniques by a small margin. However when 
the signal-noise ratio was decreased (-6 dB), the multiresolution decomposition achieved 
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superior results. Among the multiresolution structures, the results varied with the number of 
subchannels implemented and the region of the spectrum considered. 

Finally, in examining the application of multiresolution decomposition to dynamic signals, the 
results were entirely dependent on the region of the frequency spectrum in which the process 
being detected occurred. In the higher regions of the spectrum, the ability of the multiresolution 
techniques to resolve proximate narrowband components was poor at best. In the lower end of 
the spectrum, however, good results were achieved. 

B. RECOMMENDATIONS FOR FURTHER STUDY 

The most obvious direction in which to continue the study of multiresolution decomposition 
techniques is to improve the quality of the filter banks used for subchannel decomposition. In 
the structures implemented, decomposition into more than two subchannels was accomplished 
with pseudo-QMF banks. As indicated, methods for the design of perfect reconstruction filter 
banks of an arbitrary number of channels exist. The question arises as to whether results might 
not improve through the use of perfect reconstruction filter banks for decomposition into 
subchannels. 

Further study directed toward the determination of the optimum decomposition structure may 
also be warranted. In section 6.C, within the context of generation of comparing receiver 
operator characteristics for different structures, it was suggested that structures involving the 
decomposition of detail sequences into multiple subchannels may result in excessive resolution 
at low frequencies. As an alternative to decomposing the detail sequence into subchannels at 
each stage, limiting the performance of subchannel decomposition to the first two to three stages 
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may improve results. Alternatively, varying at each stage the number of subchannels into which 
the detail sequence is decomposed may also provide an option. For instance, a possible scheme 
would involve decomposition of the first stage detail sequence into three subchannels, the second 
stage detail sequence into two subchannels and performing no subchannel decomposition for the 
detail sequence at all subsequent stages. 

Additionally, the performance of multiresolution techniques for detecting dynamic signals in the 
presence of noise was not addressed. The periodogram decompositions of the quadratically 
chirped signals produced, through a spreading over multiple spectral bins, a significant reduction 
in the level of signal indication for the dynamic signal component. The steady-state harmonic 
component of the test sequence remained confined to one or two spectral bins while the dynamic 
component appeared spread over 12-16 spectral bins. Consequently, the signal levels of 
equal-power dynamic and steady-state components of the signal varied by several decibels. This 
did not occur in the case of the multiresolution decomposition. Since, for the multiresolution 
decomposition the dynamic signal and the steady-state components remained equally well 
localized, the peaks of the dynamic signal component maintained the same approximate height 
as the steady-state component. It stands to reason, therefore, that the multiresolution detector 
may be less susceptible to noise when applied to the problem of a dynamic signal. 

Finally, further work is merited in the determination of the optimum manner in which to present 
multiresolution decompositions. For, a given signal, the essential information provided by 
multiresolution techniques is contained in a non-uniform lattice of nodes representing signal 
components spread, in varying degrees, over time. The method demonstrated in the present 
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work was to decompose the channel into channels and subchannels and then re-expand each 
channel individually in such a manner that the sum of the contents of each channel provides the 
reconstructed signal. This process of presentation introduces a significant computational burden. 
As was indicated in the section 6.B, the number of multiplications necessary to re-expand the 
signal always exceeded the number of computations for the actual decomposition by a 
significant margin. Brooks [33], on the other had presented only the coefficients of expansion in 
a surface plot with varying numbers of zeros inserted between nodes of information at different 
scales. This approach resulted in plots such as were presented in the present study for the a trous 
algorithm. Not only does this approach provide plots from which information is relatively 
difficult to extract, but as addressed in chapter 4, the issue of delays which are disparate among 
channels arises. Further consideration regarding the best manner in which to present the results 
of multiresolution decomposition is, therefore, merited. 
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APPENDIX A DETAILS OF DEVELOPMENT OF MALLAT'S 

ALGORITHM 

In chapter 4, it was demonstrated that, since the scaling function vector lies within the 

span of the set of scaling function vectors, n } n € z , any vector <J) m+1< n can be represented as a 
Fourier series expansion in 0 m n : 

0m+l ,n(0 = T ^Qm+1. ni 0 m, k^0m. k(t)- (A.l) 

Now the series coefficient, the inner product term contained in the summation (A.l), equals, by 
definition 

(<t>m+l.n, 0m.k j = J 0m+l,n(t) ’ 0m,k(O dt. (A.2) 

Substituting (4.14) into the right-hand side of (A.2) produces 

S 0m + i.„(t)-0m.k(t)dt = 2 _(m+1)/2 2~ m2 j 0(2-< m+1 >-t-n)-0(2- m -t-k)dt 

= -4-2" m } 0(2^ m+1) -t-n) 0(2' m 't-k)dt (A3) 

= -L- 2~ m j 0( 2 ~"*- 2 '" ) • 0(2~ m • t - k) dt 

*2 

Now, (A.l) will become much simpler if the series coefficients can be expressed independent of 
the scale index m. Application of the change of variable of integration u = 1 m -t - 2 n to the 
final line of (A. 3) produces this result: 

j 0m +1 .n(t)-0 m .k(t)dt = -\=r 2~ m j 0(2^4=12.) • 0(2" m • t - k) dt 

v *- 

= j 0^yj-0^(u + 2-n)-kj du . (A.4) 

= -j I 0(f)'0(u-(k-2-n)j du 
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Since, therefore, by (4. 14), 9m.n(t)=2 m 2 9(2 m • t-n), (A. 4) becomes 



-=* ! $(7 J-9^u-(k-2-n)Jdu = -=- / <|>i.o(u)-<|)o.k-2n(u) du 

= ~=* 9 l.U> 9 <).k- 2 n 



(A. 5) 



Consequently, (4.19) is proved. From (A.5) and (4.19), (4.20) follows: 



9 m+l.n(t) — — £ ( 9 1.0, 9 1 ). k-2 n ) 9 m.k(t)- 



Developing (4.23) from (4.20) involves, first, substitution of the definition (4.14): 



2 -(m+l) ' 2 9(2" (m+u • t - n) = 2 -,n2 -=r S ^,. 0> <]) 0 , k _ 2 J 9(2“ m • t - k). 

Now, contained on each side of (A. 6 ) are factors of 2‘ m and 2’ 1 2 . Cancellation of those 



(A. 6 ) 



common factors produces 



9(2-^-t-n) = I( 9 , 0 , 4 >o.k- 2 n ) 9(2" m -t-k). 



(A. 7) 



Next, again introducing the change of variables u = 2" m t - 2n results in 

<t>(y j = S ^9 1.0, 9 o.k -2 n j 9 ^u - (k - 2 • n) j . 
Finally, translation of the index of summation k = k' - 2-n in (A. 8 ) yields 



(A. 8 ) 



9 ( 7 ] =S [ 9i. 0) 9o.k ] 9 (u-k). 



(A. 9) 



To generate the FIR filter defined in (4.23), it must first be specified that the scaling functions 9 
are compactly supported. Then, integrating both sides of (A. 9) over their domains results in 



J 9(7 ) du = I (9i.o, 9o.k ) j 9(u-k) du. (A. 10) 
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Application of the changes in variables of integration v = u/2 and vv = u - k to (A. 10) produces 



2 • j 0 (v) dv = I 01 .,;, (pi*. k j <p(w) dvv 



(A. 1 1 ) 



If the common factors of j <p( v) dv are removed from each side of (A. 1 1), the result becomes 



2 - X [ <J>i.o, 0'i.k ) • 



(A. 12) 



If h(k)=f 



<?> 0 ,k 



, then 



Sh(k)= 1 . 

k 
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APPENDIX B MA TLAB SOURCE CODE 



The source \ fat lab code for primary algorithms employed during this study are provided in this 
appendix. 



A. LAPLACIAN PYRAMID 

% 

% This MATLAB function performs FULL-CHANNEL DECOMPOSITION 
% OF A SEQUENCE IN ACCORDANCE WITH THE LAPLACIAN PYRAMID. 
% The Laplacian Pyramid algorithm is described in Michael 
% Unser, "An improved least squares Laplacian pyramid ... ", 

% Signal Processing 27 ( 1992) 187-203. 

% 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



The function syntax is: 

P=Iappmd(s, N) 

where: "s" is a vector containing the sequence 
to be decomposed. 

"N" is an integer indicating the number 
of levels of decomposition to be 
performed. 

"P” is an array whose rows contain the 
sequence component corresponding 
to that channel. 

The array "P" is computed such that summing along each 
column produces the reconstructed sequence. 



function P=lappmd(s, N) 

% 



% The decomposition kernel is defined as "w." 

0 / 

/o 



w=[l/4-.2 1/4 .4 1/4 1/4-.2]; 

% 

% The input sequence "s" is reshaped as a row vector. 

% 



s=reshape( s, 1 Jength(s)); 

% 

% The content is "s" is assigned to vector "pkl." The array 
% "D" is initialized as a zero vector. 



% 



pkl=s; 

% 



% Decomposition of the sequence is performed. 
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% 

for k=l:N 

% 



% 



% 



% 



% 

% 

% 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



% 

end 

% 

% 

0 / 

/o 

0 / 

/O 



"pic" is reinitialized for next iteration. 
pk=pkl : 

The sequence is reduced using user-defined function REDUCE. 
pkl=reduce(pk,w); 

"pkl" is symmetrically windowed to one-half the length of "pk." 
pkl=pkl(2:length(pkl)-l); 

The sequence "pkl" is expanded and the result assigned to "dk." 
dk=expand(pkl,2*w); 

"dk" is symmetrically windowed to the length of "pk." 
dk=dk(3:length(dk)-2); 

The sequence "dk" is assigned to the kth row of array "D." The 
length of "dk" is stored in the kth element of "L." 

L(k)=min(length(dk),length(pk)); 

D( k, 1 :L(k))=pk( l :L(k))-dk( 1 :L(k)); 



The final approximation is assigned to the "k+l'th" row 
of "D." 



L( length! L)+ 1 )=length( pk 1 ); 
D(length(L),l:L(length(L)))=pkl; 



0 / 

/o 

% The full-channel expansion is calculated from the 

% rows of array "D" and the results assigned to the rows 
% of "P." The first row of "P" is equal to the first row 
% of "D." 

% 

P( 1,:)=D(1,:); 

% 

for k=2:N+l 
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% 

% The detail from channel "k" is assigned to vector "dk." 

% 

dk=D(k. l:L(k)); 

% 

% The channel content is expanded. 

0 / 

/ 0 

for j=2:k 

dk=expand(dk,2*w); 
dk=dk(3: length! dk)-2); 

end 

% 

% The expanded channel content is assigned to row "k" of 
% array "P." 

% 

P(k, 1 : length! dk))=dk; 

0 / 

/() 

end 

B. A TROUS ALGORITHM 

% 

% This MATLAB function calculates an approximation of the 
% DISCRETE WAVELET TRANSFORM of a function 
% using SHENSA’S ATROUS ALGORITHM discribed in "The 

% Discrete Wavelet Transform: Wedding the A Trous and 

% Mallat Algorithms," IEEE Transactions on Signal 
% Processing, Vol. 40, No. 10, Oct. 1992. 

% 

% The function syntax is: 

% W=shendwt(s, M, P_i, P_f) 

% where: "s" represents a vector of data to be transformed 
% "M" is an integer indicating the number of "voices" 

% tu be used to cover the frequency spectrum 

% "W" is an array containing the coefficients of 

% the magnitude-squared wavelet transform 

% of”s." 

% "P_i" represents the first wavelet transform scale 

% to be calculated 

% "P_f represents the final wavelet transform scale 

% The transform is calculated using a madulated Gaussian 
% window as an analyzing wavelet. If multiple voices are 

% specified, the projection of "s" on each voice will be 

% represented separately. The Lagrange, a trous interpolation 
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% filter used is obtained from convolving a four-point 
% Daubechies scaling filter with itself, 
function W=shendwt(s, M, P_i, P_f) 

% 

% First, the argument vector "s" is conditioned. If"s"is 
% defined as a column vector, it is converted to a row vector. 

% Secondly, ”s" is zeropadded to the next integer power of "2." 

% 

[rows,cols]=size(s); 

% 

if cols == 1 
s=s.'; 

end 

clear rows;clear cols; 

0 / 

/O 

s=zeropad(s, l,2 A ceil(log(length(s))/log(2))); 

% 

% Default values are imposed for starting and ending scales if 
% none are specified. 

% 

if exist('P_i') == 0 

PJ=l; 

end 

if exist('P_f) == 0 

P_f=floor(log(length(s))/log(2)); 

end 

% 

% If the number of voices is not specified, a default value of 
% M=2 is imposed. 

% 

if exist('M') == 0 
M=2; 

end 

0 / 

/() 

% Next the analyzing wavelet must be calculated. The starting 
% point of this process is to specify the Gaussian window rolloff 
% factor "beta" in accordance with the specified number of 
% voices. (Shensa (6.31)). If M=l, the value of "beta" is defined 
% as "pi/(4*sqrt(2))." 

% 

if M == 1 

% 

beta=pi/(4*sqrt(2)); 

% 
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else 

% 

beta=l/(2*M); 

% 

end 

% 

% The location of the center frequency "nu" is assigned in accordance 
% with Shensa (6.27). 

% 

nu=pi-sqrt(2)*beta; 

% 

% The voice scaling factor "a" is calculated using Shensa (6.32). 

% 

a=2 A (l/M); 

% 

% The region of support for the analyzing wavelet filter "g” is 
% approximated as the region for which the Gaussian window 
% is greater than 10 A -3. Consequently, the filter impulse 
% response domain is [-a A (M-l)*sqrt(14)/beta, a A (M-l)*sqrt( 14)/beta], 
% The length of "g" is represented by an odd integer. 

% 

n=-ceil(a A (M-l)*sqrt(14)/beta):ceil(a A (M-l)*sqrt(14)/beta); 

% 

% The analyzing wavelet is calculated for each voice in accordance 
% with Shensa (6.32). It is then normalized such that its peak 
% passband frequency response is unity. 

% 

for k=l:M 
% 

g(k,:)=exp(j *nu*(n/(a A (k- 1 )))). *exp(-(beta A 2 *(n/a A (k- 1 )). A 2)/2); 

% 

end 

clear n; 

% 

% A variable "G" indicating the half-length of filter "g" is 
% defined. 

% 

G=ceil(length(g)/2); 

% 

% Next, the Lagrange interpolation filter is calculated. The filter 

% is obtained from "auto-convolution" of a Daubechie four-point 

% DWT filter. 

% 

f=[l+sqrt(3) 3+sqrt(3) 3-sqrt(3) l-sqrt(3)]/(4*sqrt(2)); 
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% 

f=conv(f,fliplr( f)); 
f=f/sqrt(2); 

% 

% A variable "F" indicating the half-length of "P is defined. 

0 / 

' 0 

F=ceil(length( f), 2); 

% 

% The recursion is next executed according to Shensa (2. 12a & b). 
% The sequence is filtered and decimated the first "P_i" times. 

% 

for k=l:P_i-l 

% 

s=conv(s,f); 

% 

s=s(F:2:length(s)-F); 

% 

end 

% 

% The output matrix "W" is initialized as a zero vector of 
% dimensions identical to those of "s." 



% 

W=zeros(s); 

% 



for k=P_i:P_f 

% 

% The data vector "s" is first filtered with each voice of "g." 

% The squared magnitude of the result is assigned to the appropriate 

% column of "W." 

% 

for n=(k-P_i)*M+l:(k-P_i+l)*M 

% 



% The row of "W" to be evaluated is initialized as zero. 

% 

W(n,:)=zeros(W(l,:)); 

% 

% The squared magnitude of the filter output is calculated 
% and assigned to "Wk." 

% 

Wk=sqrt(abs(conv(s,g(n-(k-P_i)*M,:))). A 2); 

% 

% The elements of Wk" are assigned to the corresponding 
% columns of the row of "W" being calculated. 

% 
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% 

% 

% 



% 

% 

end 



W(n,2 A (k-P_i):2 A (k-P_i):length(W))=Wk(G:length(Wk)-G+l); 

end 

"s" is filtered with the lagrange interpolation filter and then 
decimated. 

s=conv(s,f); 
s=s(F:2: length(s)-F); 



C. MULTIRESOLUTION DECOMPOSITIONS FROM CASCADES OF FILTER 
BANKS 

% 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



/ 0 

% 

% 

% 

% 

% 

% 



This MATLAB function performs WAVELET-LIKE DECOMPOSITION 

of a sequence using CASCADES OF QMF-TYPE FILTER BANKS 

and then performs FULL EXPANSION BY CHANNEL 

for each channel of decomposition. At each stage, the 

sequence is decomposed into an approximation channel and 

one or more detail channels. The detail channels can then 

be further decomposed into subchannels. 

The function syntax is: 

D=wavebank(s,h,N,g) 

where: "s" is a vector representing the sequence to 
be decomposed 

"h" is an array whose rows contain the filter 
coefficients for the primary analysis filter 
bank in order of increasing center frequency. 

"N" is an integer indicating the number of stages 
of decomposition to perform. 

"g" is an array whose rows contain the filter 
coefficients for the analysis filter for the 
subchannel decomposition in order of 
increasing center frequency. (OPTIONAL) 

If "g" is not specified, the function, by default, does not 
decompose the detail channel into subchannels. 
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function D=wavebank(s,h,N,g) 

% 

% The dimensions of "h" and "g" (if defined) are stored in 
% variables Rc, Cc, Rs, Cs, indicating the array dimensions for 
% channel and subchannel decompositions, respectively. The number 
% of rows will be used as the decimation and expansion factors. 

% Then, the analysis filter banks are defined and stored in 
% "f for the primary channels and "i" for the subchannels. 

% Because of the effects of decimation, the synthesis filters 
% must be scaled by a the decimation factor. 

% 

[Rc,Cc]=size(h); 

f=Rc*fliplr(h); 

Lc=length(h)-1; 

% 

Ls=0; 

Rs= 1 ; 

if exist('g’) == 1 

% 

[Rs, Cs]=size(g); 
i=Rs*fliplr(g); 

Ls=length(g)-1; 

% 

end %if 

% 

% The sequence "s" is reshaped as a row vector. 

% 

s=reshape(s, 1 ,length( s)); 

% 

% The sequence is decomposed. At each stage, the 
% approximation of "s" is stored in a vector "sk." 

% The next stage of decomposition is stored in "ski." 

% 

$kl(Rc,:)=s; 

% 

for k=0:N-l 

% 

% Sequence vector s is reinitialized for next iteration. 

% 

sk=skl; 

() / 

/O 

% The sequence "s" is decomposed into primary channels by 
% the function REDUCE and the result is assigned to "ski." 

% 
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for p=l:Rc 

% 

dp=reduce( sk( Rc,: ),h( Rc-p+ 1 ),Rc); 
skl(p, 1 :length(dp))=dp; 

% 

end %p 

% 

% "ski" is truncated to its proper post-"REDUCE” length. 

% 

sk 1 — sk l ( :, 1 : length(dp)); 

% 

% If "g" is defined, the sequence "dkl" is reduced into Rc 
% subchannels. The lengths of the subchannel sequences for 
% each channel are stored in a vector "L." 

% 

for p=l:Rc-l 
if exist('g’) == 1 

% 

for q=l:Rs 

% 

dq=reduce(skl(p,:),g(Rs-q+l,:),Rs); 

L(k+l)=length(dq); 

D((Rc-l )*Rs*k+Rs*(p- l)+q, 1 :L(k+ 1 ))=dq; 

% 

end %q 

% 

else 

% 

L(k+l)=length(dp); 

D((Rc-l)*k+p,l :L(k+l))=skl(p,:); 

% 

end %if 
end %p 

% 

end %k 

% 

keyboard 

% 

% The sequence is reconstructed using full expansion by channel. 

% 

for k=l:N 
% 

% If "g" is defined, the detail "d" for each subchannel is 
% extracted from the corresponding row of "D." 
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% 



if exist(’g') == 1 

% 

for p= 1 : Rc- 1 
for q=l:Rs 

% 

d=D(Rs*( Rc-1 )*(k-l)+Rs*(p-l)+q, 1 :L( k)); 

% 

% The detail is expanded with the corresponding subchannel 
% synthesis filter. 

% 

d=expand(d,i(Rs-q+l,:),Rs); 

% 

% The subchannel component is expanded with the primary detail 
% channel synthesis filter "f(@,:)" and then by the approximation 
% channel synthesis filter "f( 1,:)" the number of times 
% appropriate for the number of stages of decomposition performed. 

% 

d=expand(d,f(Rc+l-p,:),Rc); 

% 

for m=l:k-l 



% 

% 

% 

% 

% 

% 

% 

% 



% 



% 

% 

% 

% 

% 

else 



d=expand(d,f( l,:),Rc); 



end %m 

The total delay resulting from transmission through 
the system is calculated and assigned to "delay." It 
is then removed and the expanded sequence "d" is 
reintroduced to the corresponding row of array "D." 

delay=sum(fRc). A [0:k-l])*Lc+Rc A k*Ls; 

d=d(delay+l:delay+length(s)); 

M(k)=length(d); 

D(Rs*(Rc- 1 ) *(k- 1 )+Rs*(p- 1 )+q, 1 : M(k))=d; 

end %q 
end %p 

If "g" is not defined, the subchannel reconstruction steps 
are bypassed. Reconstruction of the primary detail channel 
is performed. 
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% 



% 

% 

% 

% 

% 

% 



% 



% 



% 



for p=l:Rc-l 

d=D( ( Rc- 1 )*(k- 1 )+p, 1 : L(k)); 

The detail sequence is expanded by the primary detail channel 
synthesis filter "f(Rc+l-p,:)" and then by the approximation 
channel synthesis filter "f( 1 )" as many times as 
appropriate for the number of stages of decomposition performed. 

d=expand(d,f(Rc+l-p,:),Rc); 

for m=l:k-l 

d=expand(d,f( 1 ,0,Rc); 

end %m 



% 

% The total delay resulting from transmission through 
% the system is calculated and assigned to "delay." It 
% is then removed and the expanded sequence "d" is 
% reintroduced to the corresponding row of array "D." 

% 

delay=sum((Rc). A [0:k-l])*Lc; 
d=d(delay+l :delay+length(s)); 

M(k)=length(d); 

D((Rc- 1 )*(k- 1 )+p, 1 :M(k))=d; 

% 

end %p 

% 

end %if 



% 

end %k 



0 / 

/O 

% The final approximation sequence is selected as the final 
% row of "ski." 



% 

skl=skl(Rc,:); 

% 

% The final approximation signal is expanded with primary 
% approximation channel synthesis filter according to the 
% number of decomposition stages. 

% 

for m=l:N 

% 
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skl=expand(skl,f(l,:),Rc); 



% 

end %m 

% 

% The expanded approximation sequence is added as the final row 
% of array "D." 

% 

% The total delay resulting from transmission through 
% the system is calculated and assigned to "delay." It 
% is then compensated and the expanded sequence "ski" is 
% reintroduced to the final row of array "D." 

% 

delay=sum((Rc). A [0:N-l])*Lc; 

skl=skl(delay+l:delay+length(s)); 

M( length! M)+ 1 )=length(sk 1 ); 

D( Rs*(Rc- 1 )*N+ 1 , 1 :M(length(M)))=sk 1 ; 



D. PERIODOGRAM DECOMPOSITION 

% 

% This MATLAB function calculates the PERIODOGRAM SPECTRAL ESTIMATE 
% for a specified input data vector. The function syntax is 

% 

% S=pgram(s,w, N); 

% 

% where: "s" is the data vector to be analyzed. 

% "w" represents a window with which data will be windowed 

% "N" is an integer indicating the number of steps to 

% shift the window for each step 

% (default = length(w)/2) 

% "S" is a periodogram-type output matrix. 

0 / 

/o 

% The argument data vector is subdivided with 50% overlap. Each 
% row of the output array represents the transform of a subdivision 
% of the original data vector. 

% 

function S=pgram(s,w,N); 

% 

% Check to see if "w" is defined. If "w" is not defined, 

% a rectangular window of length "length(s)" is imposed as 
% a default. 

% 

if exist('w’) == 0 
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% 

w=ones(s); 

% 

end 

% 

% A default vale is imposed for "N" 

% 

if exist('N') == 0 

% 

N=length(w)/2; 

% 

end 

% 

% Convert column vector data vector arguments to row vector 
% format. 

% 

s=reshape(s,length(s),l); 
w=reshape( w,length( w), 1 ) ; 

% 

% "s" is zeropadded to an integer number of window lengths. 

% 

for k= 1 :(length(s)-length(w))/N+ 1 

% 

r=fft(w.*s((k- 1 )*N+l:(k- l)*N+length(w))); 

% 

S(:,k)=r( 1: length! w)/2); 

% 

end 

0 / 

/o 

% In accordance with the definition of the periodogram, the 
% periodogram magnitude is divided by the window length. 

% 

S=S/length(w); 

E. DYADIC EXPANSION OF TWO-SCALE DIFFERENCE EQUATIONS 

% 

% This function expands a two-scale, difference equation of 

% the form 

% phi(x)=sum(c_k*phi(2*x-k)). 

% The function syntax is 

% [phi, psi]=twoscale(c,P,eigtol) 

% where the function arguements are: 

% "c" is a coefficient vector [c_0 c_l ... c_N-l] 

% "P" is the number of dyadic points evaluated between 
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each integer in the domain. 

"eigtol" is indicates the maximum deviation from unity 
which will be tolerated in calculating the eigenvalue 
corresponding to the eigenvector containing the 
values of "phi" for integer arguments. 

(DEFAULT = IE-3) 

The output variables are 

"phi" represents the scaling function satisfying the 
difference equation 
"psi" is the corresponding wavelet. 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

function [phi,psi]=twoscale(c,P, eigtol) 

% 

% If "eigtol" is not defined, a default vale is imposed. 

% 

if exist(’eigtol') ~= 1 
eigtol=le-3; 
end 

% The argument vector "c" is reshaped as a column vector. 

% 

c=reshape(c,length(c), 1); 

% 

% The condition that the sum of the elements of "c" equals "2" 
% is imposed. 

% 

c=2*c/sum(c); 

0 / 

/o 

% First the scale function is evaluated for integer points 
% on its domain. A zero-padded version "c_n" of the coefficent 
% vector "c" is created. 

% 

c_n=[zeros(l:length(c)-l)'; c; zeros(l:length(c)-l)']; 

% 

% A matrix "C" containing values of the coefficient vector 
% "c" is generated. 



for k=l:length(c) 

% 

C(k,:)=flipud(c_n(2*k- 1 :2*(k- 1 )+length(c)))'; 

% 

end 

% 

% The matrix "C" is truncated to eliminate rows corresponding 
% to the trivial solutions corresponding to non-unity 
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% "c_0" and "c_N-l". 

% 

if c(l) ~= 1 

(> /o 

C=C(2:length(C),2:length(C)); 

% 

end 

% 

if c(length(c)) 1 

% 

C=C( 1 : length(C)- 1 , 1 : length(C)- 1 ); 

% 

end 

% 

% The resulting matrix is checked to ensure that it contains 
% at least one eigenvalue of unity. If one unity eigenvalue 

% is present, the function continues. Otherwise, it terminates 

% in an error message. 

% 

% The eigenvectors and eigenvalues of "C" are calculated. 

% The diagonal matrix containing the eigenvalues is converted to 
% a column vector of eigenvalues. 

% 

[E,D]=eig(C); 

D=diag(D); 

% 

if min(abs(D-ones(D))) <= eigtol 

% 

% The address of the unit eigenvalue is determine. The 
% eigenector corresponding to that address is defined as "phi." 

% Values of "zero" are imposed on the endpoints for non-unity 
% "c_0" or "c_N-l". 

% 

K=find(abs(D-ones(D)) == min(abs(D-ones(D)))); 

K=K(1); 

phi=E(:,K); 

% 

ifc(l) ~= 1 

% 

phi=[0:phi]; 

% 

end 

% 

if c(length(c)) ~= 1 
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% 



phi=[phi;0]; 



% 



% 

% 

% 

% 

% 

% 

% 



% 

% 

% 

% 

% 

% 

% 

% 

% 



end 



The values of the scaling function are caluclated recursively 
for dyadic values of "phi." The values of "phi" at half-integer 
domain points is merely the convolution of "phi" and "c." 

N=length(c)-l; 
for k=0:max(0,P-l) 

A circulant, convolution matrix "M" is created. 

The matrix "M" is cleared from memory for the next iteration, 
clear M 
for m=0:N 

M(:,m+l)=[zeros(m*2 A k, l);phi;zeros((N-m)*2 A k, 1)]; 
end 



%size(M) 

% The convolution is evaluated. 

% 

phi=M*c; 

% 



% 



end 

% 

% The corresponding wavelet is calculated. 

% 

n=0:length(c)-l; 

c=flipud(c).*((-l). A n'); 

psi=M*c; 

% 

% An error message is displayed in the event no 
% eigenvalues of unity are present. 

% 

else 
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% 

'Filter string produces no unity-value eigenvalues.' 

% 

end 



H. GENERATION OF RECEIVER OPERATING CHARACTERISTICS 

% This MATLAB script file generates receiver operating characteristics 
% for various multiresolution decompositions, saves the results 
% to files, and when complete, terminates MATLAB. 

% 

% Generate roc curve for two-channel, four-subchannel, 

% SNR = - 3 dB, high-frequency transient 

% 

load slf 
load pqmfb4 
g4=h; 

load wavebank 
D=wbdec(s,h,8,g4); 

K=find(D. A 2 >= max(max(D. A 2))/100); 

rand('normal’); 

sig=sqrt(mean(s. A 2)); 

% 

'two-channel, four-subchannel' 
for k= 1:500 

n=4*sig*rand(s); 

D=wbdec(n,h,8,g4); 

H0(k)=sum(sum(D(K). A 2)); 

D=wbdec(n+s,h,8,g4); 

H 1 (k)=sum(sum(D(K). A 2)); 
end 
% 

[PFA231f,PD231f]=roc(H0, HI, 50); 
save roc231f PFA231f PD231f 

% 

clear 

% 

clear 

% 

% Generate roc curve for two-channel, three-subchannel, 

% SNR = - 3 dB, high-frequency transient 

% 
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load slf 

load wavebank 

D=\vbdec(s,h,8,g); 

K=find(D. A 2 >= max(max(D. A 2))/100); 

rand('normal'); 

sig=sqrt(mean( s. A 2 ) ): 

% 

'two-channel, three-subchannel' 
for k= 1:500 

n=4*sig*rand(s); 

D=wbdec(n,h,8,g); 

H0(k)=sum(sum(D( K). A 2)); 

D=wbdec( n-f-s,h,8,g); 

Hl(k)=sum(sum(D(K). A 2)); 

end 

% 

[PFA231f,PD231f]=roc(H0, HI, 50); 
save roc231f PFA231f PD231f 

% 

clear 

% 

% Generate roc curve for two-channel, two-subchannel, 
% SNR = - 3 dB, high-frequency transient 

% 

load slf 

load wavebank 

D=wbdec(s,h,8,h); 

K=fmd(D. A 2 >= max(max(D. A 2))/100); 
randCnormal’); 

sig=sqrt(mean(s. A 2)); 

0 / 

/O 

'two-channel, two-subchannel' 
for k= 1:500 

n=4*sig*rand(s); 

D=wbdec(n,h,8,h); 

H0(k)=sum(sum(D(K). A 2)); 

D=wbdec(n+s,h,8,h); 

H 1 (k)=sum(sum(D(K). A 2)); 
end 
% 

[PFA221f,PD221f]=roc(H0, HI, 50); 
save roc221f PFA221f PD221f 

% 
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clear 

% 

% 

% Generate roc curve for two-channel, zero-subchannel, 
% SNR = - 3 dB, high-frequency transient 

% 

load slf 

load wavebank 

D=wbdec(s,h,8); 

K=find(D. A 2 >= max(max(D. A 2))/100); 
randf normal'); 
sig=sqrt(mean(s. A 2)); 

% 

'two-channel, zero-subchannel' 
for k= 1:500 

n=4*sig*rand(s); 

D=wbdec(n,h,8); 

H0(k)=sum(sum(D(K) A 2)); 

D=wbdec(n+s,h,8); 

Hl(k)=sum(sum(D(K.). A 2)); 

end 

% 

[PFA201f,PD201f]=roc(H0, HI, 50); 
save roc201f PFA201f PD201f 

% 

clear 

% 

% 

% Generate roc curve for three-channel, zero-subchannel, 
% SNR = - 3 dB, high-frequency transient 

% 

load slf 

load wavebank 
D=wbdec(s,g,8); 

K.=ftnd(D. A 2 >= max(max(D. A 2))/100); 

randCnormal'); 

sig=sqrt(mean(s. A 2)); 

% 

'three-channel, zero-subchannel' 
for k= 1:500 

n=4*sig*rand(s); 

D=wbdec(n,g,8); 

H0(k)=sum(sum(D(K) A 2)); 

D=wbdec(n+s,g,8); 
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H 1 ( k)=sum( sum( D( K). A 2)); 



end 

% 

[PFA30lf,PD30lf]=roc(H0, HI, 50); 
save roc301f PFA301f PD301f 

% 

clear 

% 

clear 

% 

% 

% Generate roc curve for two-channel, four-subchannel, 
% SNR = - 3 dB, high-frequency transient 

% 

load shf 
load pqmfb4 
g4=h; 

load wavebank 
D=wbdec(s,h,8,g4); 

K=find(D. A 2 >= max(max(D A 2))/100); 

rand('normar); 

sig=sqrt(mean(s. A 2)); 

% 

'two-channel, four-subchannel' 
for k=l :500 

n=4*sig*rand(s); 

D=wbdec(n,h,8,g4); 

H0(k)=sum( sum( D( K). A 2)); 

D=wbdec(n+s,h,8,g4); 

Hl(k)=sum(sum(D(K). A 2)); 

end 

% 

[PFA23hf,PD23hf]=roc(H0, HI, 50); 
save roc23hf PFA23hf PD23hf 

% 

clear 

% Generate roc curve for two-channel, two-subchannel, 
% SNR = - 3 dB, high-frequency transient 

% 

load shf 

load wavebank 

D=wbdec(s,h,8,h); 

K=find(D. A 2 >= max(rnax(D. A 2))/100); 
randCnormal'); 
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sig=sqrt(mean( s. A 2)); 

% 

'two-channel, two-subchannel' 
for k= 1:500 

n=4*sig*rand(s); 

D=wbdec(n,h.8,h); 

H0( k)=sum( sum( D(K). A 2)); 

D=wbdec( n+s,h,8,h); 

H 1 (k)=sum( sum(Df K). A 2)); 

end 

% 

[PFA22hf,PD22hf]=roc(H0, HI, 50); 
save roc22hf PFA22hf PD22hf 

% 

clear 

% 

% 

% Generate roc curve for two-channel, zerp-subchannel, 
% SNR = - 3 dB, high-frequency transient 

% 

load shf 

load wavebank 

D=wbdec(s,h,8); 

K=find(D. A 2 >= max(max(D. A 2))/100); 

rand('normar); 

sig=sqrt(mean(s A 2)); 

% 

'two-channel, zero-subchannel' 
for k= 1:500 

n=4*sig*rand(s); 

D=wbdec(n,h,8); 

H0(k)=sum(sum(D(K) A 2)); 

D=wbdec(n+s,h,8); 

Hl(k)=sum(sum(D(K). A 2)); 

end 

% 

[PFA20hf,PD20hf]=roc(H0, HI, 50); 
save roc20hf PFA20hf PD20hf 
% 

% 

quit 

% 
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G. GENERAL-PURPOSE ROUTINES CALLED BY OTHER ROUTINES 

% 

% This function zero-pads an argument matrix forming 
% a new matrix of the specified size. The original matrix 

% will be in the upper left-hand corner the new matrix. 

% 

% The function syntax is: 

% 

% X=zeropad(x,Nl,N2) 

% 

% where: 

% "x" is the initial argument matrix 

% "Nl" is the total number of rows for the new matrix 

% "N2" is the total number of columns for the new matrix 

% 

function X=zeropad(x,Nl,N2) 

% 

[Ml,M2]=size(x); 

% 

X=[x zeros(Ml,N2-M2)]; 

% 

X=[X; zeros(Nl-Ml,N2)]; 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



This MATLAB function reduces A SEQUENCE IN THE SENSE OF 
THE EXPAND OPERATION EMPLOYED IN LAPLACIAN PYRAMID 
DECOMPOSITION. The operation entails FIR filtering 
followed by decimation in time. 

The function syntax is 

fkl=reduce(fk, h, M) 

where: fk is a vector containing the sequence to be expanded 

h is a vector containing the FIR filter coefficients 

M is an integer indicating the decimation 
factor (default == 2) 

Regardless of the shape of the input vectors, output vectors 
are returned as row vectors. 



% 

function fkl=reduce(fk, h, M); 



% 



% A default value of two is imposed on "M.” 

% 
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if exist('M') ~= 1 
M=2; 

end 

% 

% The sequences "fk" and "h" are reshaped as row vectors. 

% 

fk=reshape( fk, 1 , length! fk)); 
h=reshape(h, 1 Jength(h)); 

% 

% The sequence "fk" is filtered with "h" and the result is 
% assigned to "fkl." 

% 

fkl=conv(fk,h); 

% 

% The sequence "fkl" is decimated by a factor of "M." 

% 

fkl=fkl(l:M:length(fkl)); 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



This MATLAB function EXPANDS A SEQUENCE IN THE SENSE OF 
THE EXPAND OPERATION EMPLOYED IN LAPLACIAN PYRAMID 
DECOMPOSITION. The operation entails FIR filtering 
followed by decimation in time. 

The function syntax is 

fkl=expand(fk, h, M) 

where: fk is a vector containing the sequence to be expanded 

h is a vector containing the FIR filter coefficients 

M is an integer indicating the decimation 
factor (default == 2) 

Regardless of the shape of the input vectors, output vectors 
are returned as row vectors. 



% 

function fkl=expand(fk, h, M); 

% 



% A default value of two is imposed on "M." 

% 



if exist('M’) ~= 1 
M=2; 



end 



<>/ 

/o 

% The sequences "fk" and "h" are reshaped as row vectors. 

% 
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fk=reshape( fk, 1 ,length( fk)); 
h=reshape( h, 1 ,length(h)); 

% 

% Zero-insertion is performed on the input sequence "fk." 
% The result is assigned to "fkl." 

% 

fkl=zeros( l,M*length(fk)-l); 
fkl(l:M:length(fkl))=fk; 

% 

% The sequence "fk" is filtered with "h" and the result is 
% assigned to "fkl." 

% 

fkl=conv(fkl,h); 

% 



% This function performs numerical integration of a 
% vector in accordance with Simpson's Rule. The function 
% arguments are: 

% 

% "x"--the domain of the argument function 

% "y"— the range of the argument function. 

% 

% The function's syntax is: 

% 

% Y=simpson(x,y) 

% 

% A regular partition is assumed. 

% 

function Y=simpson(x,y) 

% 

% Vectors are checked for row- or column-vector format. Column 

% vectors are converted to row vectors. 

0 / 

/o 

dimens=size(x); 

N=max(dimens); 
if dimens(l) == N, 
x=x'; 

y=y'; 

end 

% 

% The partition width is determined and the argument range vector 
% is multiplied by the partition width. 

% 

y=y*((max(x)-min(x))/max(size(x))); 
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The integral is evaluated. 



% 

% 

% 

Y( 1 j=y(l); 
for k=2:N 

Y(k)=y( l:k)*ones(y( l:k)'); 

end 

% 



% 

% Given two, equal-length vectors, this MATLAB function 
% GENERATES HISTOGRAMS AND CONVERTS THEM TO RECEIVER 

% OPERATING CHARACTERISTICS. 

% 

% The function syntax is: 

% [PFA, PD]=roc(H0, HI) 

% where: "HO" is a vector containing "hypothesis false" 

% realizations of a random process. 

% "HI" is a vector containing "hypothesis true" 

% realizations of a random process. 

% "PFA" is a vector representing the probability 

% of false alarm 

% "PD" is a vector representing the probability 

% of detection. 

% 

function [PFA, PD]=roc(H0, HI, Nobins) 

% 

% The minimum and maximum bin locations are identified. 

% 

Bmin=min(min([HO;Hl])); 

Bmax=max(max([HO;Hl])); 

% 

% A vector is created representing the bin locations. 

<>/ 

/O 

if exist('Nobins') == 0 

% 

Nobins=length(H0)/25; 

% 

end 

B=Bmin:(Bmax-Bmin)/Nobins:Bmax; 

% 

% The MATLAB function "HIST" is invoked to produce histograms 
% of the realization vectors. 
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% 

[pfa,X]=hist(HO,B); 

[pd,X]=hist(Hl.B); 

% 

% The histograms are integrated using function "SIMPSON." 

% 

PFA=simpson( X.pfa); 

PD=simpson( X,pd); 

% 

% The integrals are normalized for maximum values of unity. 

% 

PFA=PFA/max(PFA); 

PD=PD/max(PD); 

% 

% The actual probability detect/probability false-alarm 

% plots are the complements of the probability distribution 
% functions. 

% 

PFA=1-PFA; 

PD=1-PD; 
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